From de370c35f9092fdad99af08b50cb97b7db9ca56d Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 9 Dec 2024 08:15:49 +0100 Subject: [PATCH 1/6] add seperate show_plot_makie --- src/callbacks_step/visualization.jl | 53 +++++++++++++++++++++++++---- 1 file changed, 46 insertions(+), 7 deletions(-) diff --git a/src/callbacks_step/visualization.jl b/src/callbacks_step/visualization.jl index 302e7e4462a..8b6b6f5ac13 100644 --- a/src/callbacks_step/visualization.jl +++ b/src/callbacks_step/visualization.jl @@ -13,6 +13,7 @@ mutable struct VisualizationCallback{SolutionVariables, VariableNames, PlotDataC show_mesh::Bool plot_data_creator::PlotDataCreator plot_creator::PlotCreator + figure plot_arguments::Dict{Symbol, Any} end @@ -98,7 +99,7 @@ function VisualizationCallback(; interval = 0, visualization_callback = VisualizationCallback(interval, solution_variables, variable_names, show_mesh, - plot_data_creator, plot_creator, + plot_data_creator, plot_creator, nothing, Dict{Symbol, Any}(plot_arguments)) # Warn users if they create a visualization callback without having loaded the Plots package @@ -110,8 +111,9 @@ function VisualizationCallback(; interval = 0, # Requires.jl only when Plots is present. # In the future, we should update/remove this warning if other plotting packages are # starting to be used. - if !(:Plots in names(@__MODULE__, all = true)) - @warn "Package `Plots` not loaded but required by `VisualizationCallback` to visualize results" + if !(:Plots in names(@__MODULE__, all = true) || :Makie in names(@__MODULE__, all = true)) + println(names(@__MODULE__, all = true)) + @warn "Package `Plots` or `Makie` required by `VisualizationCallback` to visualize results" end DiscreteCallback(visualization_callback, visualization_callback, # the first one is the condition, the second the affect! @@ -145,7 +147,7 @@ end function (visualization_callback::VisualizationCallback)(integrator) u_ode = integrator.u semi = integrator.p - @unpack plot_arguments, solution_variables, variable_names, show_mesh, plot_data_creator, plot_creator = visualization_callback + @unpack plot_arguments, solution_variables, variable_names, show_mesh, plot_data_creator, plot_creator, figure = visualization_callback # Extract plot data plot_data = plot_data_creator(u_ode, semi, solution_variables = solution_variables) @@ -158,7 +160,7 @@ function (visualization_callback::VisualizationCallback)(integrator) # Create plot plot_creator(plot_data, variable_names; show_mesh = show_mesh, plot_arguments = plot_arguments, - time = integrator.t, timestep = integrator.stats.naccept) + time = integrator.t, timestep = integrator.stats.naccept, figure = figure) # avoid re-evaluating possible FSAL stages u_modified!(integrator, false) @@ -168,7 +170,7 @@ end """ show_plot(plot_data, variable_names; show_mesh=true, plot_arguments=Dict{Symbol,Any}(), - time=nothing, timestep=nothing) + time=nothing, timestep=nothing, figure=nothing) Visualize the plot data object provided in `plot_data` and display result, plotting only the variables in `variable_names` and, optionally, the mesh (if `show_mesh` is `true`). Additionally, @@ -184,7 +186,7 @@ See also: [`VisualizationCallback`](@ref), [`save_plot`](@ref) """ function show_plot(plot_data, variable_names; show_mesh = true, plot_arguments = Dict{Symbol, Any}(), - time = nothing, timestep = nothing) + time = nothing, timestep = nothing, figure = nothing) # Gather subplots plots = [] for v in variable_names @@ -258,4 +260,41 @@ function save_plot(plot_data, variable_names; filename = joinpath("out", @sprintf("solution_%09d.png", timestep)) Plots.savefig(filename) end + +""" + show_plot_makie(plot_data, variable_names; + show_mesh=true, plot_arguments=Dict{Symbol,Any}(), + time=nothing, timestep=nothing, figure=nothing) + +Visualize the plot data object provided in `plot_data` and display result using Makie. +Only the variables in `variable_names` are plotted and, optionally, the mesh (if +`show_mesh` is `true`). Additionally, `plot_arguments` will be unpacked and passed as +keyword arguments to the `Plots.plot` command. + +This function is the default `plot_creator` argument for the [`VisualizationCallback`](@ref). +`time` and `timestep` are currently unused by this function. + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. + +See also: [`VisualizationCallback`](@ref), [`save_plot`](@ref) +""" +function show_plot_makie(plot_data, variable_names; + show_mesh = true, plot_arguments = Dict{Symbol, Any}(), + time = nothing, timestep = nothing, figure = nothing) + + if figure === nothing + @warn "Creating new figure" + figure = GLMakie.Figure() + end + + for v in 1:size(variable_names)[1] + GLMakie.volume(figure[intTo2DInt(v)...], plot_data.data[v]; + algorithm = :mip, plot_arguments...) + end + + GLMakie.display(visualization_callback.figure) + + # TODO: show_mesh +end end # @muladd From 74f6fb7c0f493ca09753fd0930a8a5ce156a273e Mon Sep 17 00:00:00 2001 From: s6nistam Date: Thu, 12 Dec 2024 10:46:01 +0100 Subject: [PATCH 2/6] Visualization Callback Makie current state --- src/Trixi.jl | 8 + src/callbacks_step/visualization.jl | 609 +++++++++++++++------------- 2 files changed, 333 insertions(+), 284 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 7d557ddde38..bfef3fc0dd3 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -314,6 +314,14 @@ function __init__() end end + @require GLMakie="e9467ef8-e4e7-5192-8a1a-b1aee30e663a" begin + using .GLMakie: GLMakie + end + + @require CairoMakie="13f3f980-e62b-5c42-98c6-ff1f3baf88f0" begin + using .CairoMakie: CairoMakie + end + @static if !isdefined(Base, :get_extension) @require Convex="f65535da-76fb-5f13-bab9-19810c17039a" begin @require ECOS="e2685f51-7e38-5353-a97d-a921fd2c8199" begin diff --git a/src/callbacks_step/visualization.jl b/src/callbacks_step/visualization.jl index 8b6b6f5ac13..3748f265eae 100644 --- a/src/callbacks_step/visualization.jl +++ b/src/callbacks_step/visualization.jl @@ -3,298 +3,339 @@ # we need to opt-in explicitly. # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. @muladd begin -#! format: noindent - -mutable struct VisualizationCallback{SolutionVariables, VariableNames, PlotDataCreator, - PlotCreator} - interval::Int - solution_variables::SolutionVariables - variable_names::VariableNames - show_mesh::Bool - plot_data_creator::PlotDataCreator - plot_creator::PlotCreator - figure - plot_arguments::Dict{Symbol, Any} -end - -function Base.show(io::IO, - cb::DiscreteCallback{Condition, Affect!}) where {Condition, - Affect! <: - VisualizationCallback - } - visualization_callback = cb.affect! - @unpack interval, plot_arguments, solution_variables, variable_names, show_mesh, plot_creator, plot_data_creator = visualization_callback - print(io, "VisualizationCallback(", - "interval=", interval, ", ", - "solution_variables=", solution_variables, ", ", - "variable_names=", variable_names, ", ", - "show_mesh=", show_mesh, ", ", - "plot_data_creator=", plot_data_creator, ", ", - "plot_creator=", plot_creator, ", ", - "plot_arguments=", plot_arguments, ")") -end - -function Base.show(io::IO, ::MIME"text/plain", - cb::DiscreteCallback{Condition, Affect!}) where {Condition, - Affect! <: - VisualizationCallback - } - if get(io, :compact, false) - show(io, cb) - else + #! format: noindent + + mutable struct VisualizationCallback{SolutionVariables, VariableNames, PlotDataCreator, + PlotCreator} + interval::Int + solution_variables::SolutionVariables + variable_names::VariableNames + show_mesh::Bool + plot_data_creator::PlotDataCreator + plot_creator::PlotCreator + figure + plot_arguments::Dict{Symbol, Any} + end + + function Base.show(io::IO, + cb::DiscreteCallback{Condition, Affect!}) where {Condition, + Affect! <: + VisualizationCallback + } visualization_callback = cb.affect! - - setup = [ - "interval" => visualization_callback.interval, - "plot arguments" => visualization_callback.plot_arguments, - "solution variables" => visualization_callback.solution_variables, - "variable names" => visualization_callback.variable_names, - "show mesh" => visualization_callback.show_mesh, - "plot creator" => visualization_callback.plot_creator, - "plot data creator" => visualization_callback.plot_data_creator - ] - summary_box(io, "VisualizationCallback", setup) + @unpack interval, plot_arguments, solution_variables, variable_names, show_mesh, plot_creator, plot_data_creator = visualization_callback + print(io, "VisualizationCallback(", + "interval=", interval, ", ", + "solution_variables=", solution_variables, ", ", + "variable_names=", variable_names, ", ", + "show_mesh=", show_mesh, ", ", + "plot_data_creator=", plot_data_creator, ", ", + "plot_creator=", plot_creator, ", ", + "plot_arguments=", plot_arguments, ")") end -end - -""" - VisualizationCallback(; interval=0, - solution_variables=cons2prim, - variable_names=[], - show_mesh=false, - plot_data_creator=PlotData2D, - plot_creator=show_plot, - plot_arguments...) - -Create a callback that visualizes results during a simulation, also known as *in-situ -visualization*. - -!!! warning "Experimental implementation" - This is an experimental feature and may change in any future releases. - -The `interval` specifies the number of time step iterations after which a new plot is generated. The -available variables to plot are configured with the `solution_variables` parameter, which acts the -same way as for the [`SaveSolutionCallback`](@ref). The variables to be actually plotted can be -selected by providing a single string or a list of strings to `variable_names`, and if `show_mesh` -is `true`, an additional plot with the mesh will be generated. - -To customize the generated figure, `plot_data_creator` allows to use different plot data types. With -`plot_creator` you can further specify an own function to visualize results, which must support the -same interface as the default implementation [`show_plot`](@ref). All remaining -keyword arguments are collected and passed as additional arguments to the plotting command. -""" -function VisualizationCallback(; interval = 0, - solution_variables = cons2prim, - variable_names = [], - show_mesh = false, - plot_data_creator = PlotData2D, - plot_creator = show_plot, - plot_arguments...) - mpi_isparallel() && error("this callback does not work in parallel yet") - - if variable_names isa String - variable_names = String[variable_names] + + function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{Condition, Affect!}) where {Condition, + Affect! <: + VisualizationCallback + } + if get(io, :compact, false) + show(io, cb) + else + visualization_callback = cb.affect! + + setup = [ + "interval" => visualization_callback.interval, + "plot arguments" => visualization_callback.plot_arguments, + "solution variables" => visualization_callback.solution_variables, + "variable names" => visualization_callback.variable_names, + "show mesh" => visualization_callback.show_mesh, + "plot creator" => visualization_callback.plot_creator, + "plot data creator" => visualization_callback.plot_data_creator + ] + summary_box(io, "VisualizationCallback", setup) + end end - - visualization_callback = VisualizationCallback(interval, - solution_variables, variable_names, - show_mesh, - plot_data_creator, plot_creator, nothing, - Dict{Symbol, Any}(plot_arguments)) - - # Warn users if they create a visualization callback without having loaded the Plots package - # - # Note: This warning is added for convenience, as Plots is the only "officially" supported - # visualization package right now. However, in general nothing prevents anyone from using - # other packages such as Makie, Gadfly etc., given that appropriate `plot_creator`s are - # passed. This is also the reason why the visualization callback is not included via - # Requires.jl only when Plots is present. - # In the future, we should update/remove this warning if other plotting packages are - # starting to be used. - if !(:Plots in names(@__MODULE__, all = true) || :Makie in names(@__MODULE__, all = true)) - println(names(@__MODULE__, all = true)) - @warn "Package `Plots` or `Makie` required by `VisualizationCallback` to visualize results" + + """ + VisualizationCallback(; interval=0, + solution_variables=cons2prim, + variable_names=[], + show_mesh=false, + plot_data_creator=PlotData2D, + plot_creator=show_plot, + plot_arguments...) + + Create a callback that visualizes results during a simulation, also known as *in-situ + visualization*. + + !!! warning "Experimental implementation" + This is an experimental feature and may change in any future releases. + + The `interval` specifies the number of time step iterations after which a new plot is generated. The + available variables to plot are configured with the `solution_variables` parameter, which acts the + same way as for the [`SaveSolutionCallback`](@ref). The variables to be actually plotted can be + selected by providing a single string or a list of strings to `variable_names`, and if `show_mesh` + is `true`, an additional plot with the mesh will be generated. + + To customize the generated figure, `plot_data_creator` allows to use different plot data types. With + `plot_creator` you can further specify an own function to visualize results, which must support the + same interface as the default implementation [`show_plot`](@ref). All remaining + keyword arguments are collected and passed as additional arguments to the plotting command. + """ + function VisualizationCallback(; interval = 0, + solution_variables = cons2prim, + variable_names = [], + show_mesh = false, + plot_data_creator = PlotData2D, + plot_creator = show_plot, + plot_arguments...) + mpi_isparallel() && error("this callback does not work in parallel yet") + + if variable_names isa String + variable_names = String[variable_names] + end + + visualization_callback = VisualizationCallback(interval, + solution_variables, variable_names, + show_mesh, + plot_data_creator, plot_creator, nothing, + Dict{Symbol, Any}(plot_arguments)) + + # Warn users if they create a visualization callback without having loaded the Plots package + # + # Note: This warning is added for convenience, as Plots is the only "officially" supported + # visualization package right now. However, in general nothing prevents anyone from using + # other packages such as Makie, Gadfly etc., given that appropriate `plot_creator`s are + # passed. This is also the reason why the visualization callback is not included via + # Requires.jl only when Plots is present. + # In the future, we should update/remove this warning if other plotting packages are + # starting to be used. + if !(:Plots in names(@__MODULE__, all = true) || :Makie in names(@__MODULE__, all = true)) + println(names(@__MODULE__, all = true)) + @warn "Package `Plots` or `Makie` required by `VisualizationCallback` to visualize results" + end + + DiscreteCallback(visualization_callback, visualization_callback, # the first one is the condition, the second the affect! + save_positions = (false, false), + initialize = initialize!) end - - DiscreteCallback(visualization_callback, visualization_callback, # the first one is the condition, the second the affect! - save_positions = (false, false), - initialize = initialize!) -end - -function initialize!(cb::DiscreteCallback{Condition, Affect!}, u, t, - integrator) where {Condition, Affect! <: VisualizationCallback} - visualization_callback = cb.affect! - - visualization_callback(integrator) - - return nothing -end - -# this method is called to determine whether the callback should be activated -function (visualization_callback::VisualizationCallback)(u, t, integrator) - @unpack interval = visualization_callback - - # With error-based step size control, some steps can be rejected. Thus, - # `integrator.iter >= integrator.stats.naccept` - # (total #steps) (#accepted steps) - # We need to check the number of accepted steps since callbacks are not - # activated after a rejected step. - return interval > 0 && (integrator.stats.naccept % interval == 0 || - isfinished(integrator)) -end - -# this method is called when the callback is activated -function (visualization_callback::VisualizationCallback)(integrator) - u_ode = integrator.u - semi = integrator.p - @unpack plot_arguments, solution_variables, variable_names, show_mesh, plot_data_creator, plot_creator, figure = visualization_callback - - # Extract plot data - plot_data = plot_data_creator(u_ode, semi, solution_variables = solution_variables) - - # If variable names were not specified, plot everything - if isempty(variable_names) - variable_names = String[keys(plot_data)...] + + function initialize!(cb::DiscreteCallback{Condition, Affect!}, u, t, + integrator) where {Condition, Affect! <: VisualizationCallback} + visualization_callback = cb.affect! + u_ode = integrator.u + semi = integrator.p + mesh, equations, solver, cache = mesh_equations_solver_cache(integrator.p) + if ndims(mesh) == 3 && visualization_callback.plot_data_creator == PlotData2D + visualization_callback.plot_data_creator = PlotData3D + end + visualization_callback(integrator) + return nothing end - - # Create plot - plot_creator(plot_data, variable_names; - show_mesh = show_mesh, plot_arguments = plot_arguments, - time = integrator.t, timestep = integrator.stats.naccept, figure = figure) - - # avoid re-evaluating possible FSAL stages - u_modified!(integrator, false) - return nothing -end - -""" - show_plot(plot_data, variable_names; - show_mesh=true, plot_arguments=Dict{Symbol,Any}(), - time=nothing, timestep=nothing, figure=nothing) - -Visualize the plot data object provided in `plot_data` and display result, plotting only the -variables in `variable_names` and, optionally, the mesh (if `show_mesh` is `true`). Additionally, -`plot_arguments` will be unpacked and passed as keyword arguments to the `Plots.plot` command. - -This function is the default `plot_creator` argument for the [`VisualizationCallback`](@ref). -`time` and `timestep` are currently unused by this function. - -!!! warning "Experimental implementation" - This is an experimental feature and may change in future releases. - -See also: [`VisualizationCallback`](@ref), [`save_plot`](@ref) -""" -function show_plot(plot_data, variable_names; - show_mesh = true, plot_arguments = Dict{Symbol, Any}(), - time = nothing, timestep = nothing, figure = nothing) - # Gather subplots - plots = [] - for v in variable_names - push!(plots, Plots.plot(plot_data[v]; plot_arguments...)) + + # this method is called to determine whether the callback should be activated + function (visualization_callback::VisualizationCallback)(u, t, integrator) + @unpack interval = visualization_callback + + # With error-based step size control, some steps can be rejected. Thus, + # `integrator.iter >= integrator.stats.naccept` + # (total #steps) (#accepted steps) + # We need to check the number of accepted steps since callbacks are not + # activated after a rejected step. + return interval > 0 && (integrator.stats.naccept % interval == 0 || + isfinished(integrator)) end - if show_mesh - push!(plots, Plots.plot(getmesh(plot_data); plot_arguments...)) + + # this method is called when the callback is activated + function (visualization_callback::VisualizationCallback)(integrator) + u_ode = integrator.u + semi = integrator.p + @unpack plot_arguments, solution_variables, variable_names, show_mesh, plot_data_creator, plot_creator, figure = visualization_callback + mesh, equations, solver, cache = mesh_equations_solver_cache(integrator.p) + n = Trixi.ndims(mesh) + + # Extract plot data + plot_data = plot_data_creator(u_ode, semi, solution_variables = solution_variables) + + # If variable names were not specified, plot everything + if isempty(variable_names) + variable_names = String[keys(plot_data)...] + end + + # Create plot + plot_creator(n, visualization_callback, plot_data, variable_names; + show_mesh = show_mesh, plot_arguments = plot_arguments, + time = integrator.t, timestep = integrator.stats.naccept, figure = figure) + + # avoid re-evaluating possible FSAL stages + u_modified!(integrator, false) + return nothing end - - # Note, for the visualization callback to work for general equation systems - # this layout construction would need to use the if-logic below. - # Currently, there is no use case for this so it is left here as a note. - # - # Determine layout - # if length(plots) <= 3 - # cols = length(plots) - # rows = 1 - # else - # cols = ceil(Int, sqrt(length(plots))) - # rows = div(length(plots), cols, RoundUp) - # end - # layout = (rows, cols) - - # Determine layout - cols = ceil(Int, sqrt(length(plots))) - rows = div(length(plots), cols, RoundUp) - layout = (rows, cols) - - # Show plot - display(Plots.plot(plots..., layout = layout)) -end - -""" - save_plot(plot_data, variable_names; - show_mesh=true, plot_arguments=Dict{Symbol,Any}(), - time=nothing, timestep=nothing) - -Visualize the plot data object provided in `plot_data` and save result as a PNG file in the `out` -directory, plotting only the variables in `variable_names` and, optionally, the mesh (if `show_mesh` -is `true`). Additionally, `plot_arguments` will be unpacked and passed as keyword arguments to the -`Plots.plot` command. - -The `timestep` is used in the filename. `time` is currently unused by this function. - -!!! warning "Experimental implementation" - This is an experimental feature and may change in future releases. - -See also: [`VisualizationCallback`](@ref), [`show_plot`](@ref) -""" -function save_plot(plot_data, variable_names; - show_mesh = true, plot_arguments = Dict{Symbol, Any}(), - time = nothing, timestep = nothing) - # Gather subplots - plots = [] - for v in variable_names - push!(plots, Plots.plot(plot_data[v]; plot_arguments...)) + + """ + show_plot(plot_data, variable_names; + show_mesh=true, plot_arguments=Dict{Symbol,Any}(), + time=nothing, timestep=nothing, figure=nothing) + + Visualize the plot data object provided in `plot_data` and display result, plotting only the + variables in `variable_names` and, optionally, the mesh (if `show_mesh` is `true`). Additionally, + `plot_arguments` will be unpacked and passed as keyword arguments to the `Plots.plot` command. + + This function is the default `plot_creator` argument for the [`VisualizationCallback`](@ref). + `time` and `timestep` are currently unused by this function. + + !!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. + + See also: [`VisualizationCallback`](@ref), [`save_plot`](@ref) + """ + function show_plot(ndims, visualization_callback, plot_data, variable_names; + show_mesh = true, plot_arguments = Dict{Symbol, Any}(), + time = nothing, timestep = nothing, figure = nothing) + # Gather subplots + plots = [] + for v in variable_names + push!(plots, Plots.plot(plot_data[v]; plot_arguments...)) + end + if show_mesh + push!(plots, Plots.plot(getmesh(plot_data); plot_arguments...)) + end + + # Note, for the visualization callback to work for general equation systems + # this layout construction would need to use the if-logic below. + # Currently, there is no use case for this so it is left here as a note. + # + # Determine layout + # if length(plots) <= 3 + # cols = length(plots) + # rows = 1 + # else + # cols = ceil(Int, sqrt(length(plots))) + # rows = div(length(plots), cols, RoundUp) + # end + # layout = (rows, cols) + + # Determine layout + cols = ceil(Int, sqrt(length(plots))) + rows = div(length(plots), cols, RoundUp) + layout = (rows, cols) + + # Show plot + display(Plots.plot(plots..., layout = layout)) end - if show_mesh - push!(plots, Plots.plot(getmesh(plot_data); plot_arguments...)) + + """ + save_plot(plot_data, variable_names; + show_mesh=true, plot_arguments=Dict{Symbol,Any}(), + time=nothing, timestep=nothing) + + Visualize the plot data object provided in `plot_data` and save result as a PNG file in the `out` + directory, plotting only the variables in `variable_names` and, optionally, the mesh (if `show_mesh` + is `true`). Additionally, `plot_arguments` will be unpacked and passed as keyword arguments to the + `Plots.plot` command. + + The `timestep` is used in the filename. `time` is currently unused by this function. + + !!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. + + See also: [`VisualizationCallback`](@ref), [`show_plot`](@ref) + """ + function save_plot(plot_data, variable_names; + show_mesh = true, plot_arguments = Dict{Symbol, Any}(), + time = nothing, timestep = nothing) + # Gather subplots + plots = [] + for v in variable_names + push!(plots, Plots.plot(plot_data[v]; plot_arguments...)) + end + if show_mesh + push!(plots, Plots.plot(getmesh(plot_data); plot_arguments...)) + end + + # Determine layout + cols = ceil(Int, sqrt(length(plots))) + rows = div(length(plots), cols, RoundUp) + layout = (rows, cols) + + # Create plot + Plots.plot(plots..., layout = layout) + + # Determine filename and save plot + filename = joinpath("out", @sprintf("solution_%09d.png", timestep)) + Plots.savefig(filename) end - - # Determine layout - cols = ceil(Int, sqrt(length(plots))) - rows = div(length(plots), cols, RoundUp) - layout = (rows, cols) - - # Create plot - Plots.plot(plots..., layout = layout) - - # Determine filename and save plot - filename = joinpath("out", @sprintf("solution_%09d.png", timestep)) - Plots.savefig(filename) -end - -""" - show_plot_makie(plot_data, variable_names; - show_mesh=true, plot_arguments=Dict{Symbol,Any}(), - time=nothing, timestep=nothing, figure=nothing) - -Visualize the plot data object provided in `plot_data` and display result using Makie. -Only the variables in `variable_names` are plotted and, optionally, the mesh (if -`show_mesh` is `true`). Additionally, `plot_arguments` will be unpacked and passed as -keyword arguments to the `Plots.plot` command. - -This function is the default `plot_creator` argument for the [`VisualizationCallback`](@ref). -`time` and `timestep` are currently unused by this function. - -!!! warning "Experimental implementation" - This is an experimental feature and may change in future releases. - -See also: [`VisualizationCallback`](@ref), [`save_plot`](@ref) -""" -function show_plot_makie(plot_data, variable_names; - show_mesh = true, plot_arguments = Dict{Symbol, Any}(), - time = nothing, timestep = nothing, figure = nothing) - - if figure === nothing - @warn "Creating new figure" - figure = GLMakie.Figure() + + #converts a single int into a tuple of ints, to get a square arrangement for example f(1) = (1,1) f(2) = (2,1) f(3) = (2,2) f(4) = (1,2) + function makieLayoutHelper(n) + if n == 1 + return (1,1) + end + t = intTo2DInt(n-1) + if t[1] == 1 + return (t[2] + 1, 1) + elseif t[1] > t[2] + return (t[1], t[2] + 1) + elseif t[2] >= t[1] + return (t[1] - 1, t[2]) + end end - - for v in 1:size(variable_names)[1] - GLMakie.volume(figure[intTo2DInt(v)...], plot_data.data[v]; - algorithm = :mip, plot_arguments...) + + + """ + show_plot_makie(plot_data, variable_names; + show_mesh=true, plot_arguments=Dict{Symbol,Any}(), + time=nothing, timestep=nothing, figure=nothing) + + Visualize the plot data object provided in `plot_data` and display result using Makie. + Only the variables in `variable_names` are plotted and, optionally, the mesh (if + `show_mesh` is `true`). Additionally, `plot_arguments` will be unpacked and passed as + keyword arguments to the `Plots.plot` command. + + This function is the default `plot_creator` argument for the [`VisualizationCallback`](@ref). + `time` and `timestep` are currently unused by this function. + + !!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. + + See also: [`VisualizationCallback`](@ref), [`save_plot`](@ref) + """ + function show_plot_makie(ndims, visualization_callback, plot_data, variable_names; + show_mesh = true, plot_arguments = Dict{Symbol, Any}(), + time = nothing, timestep = nothing, figure = nothing) + + if figure === nothing + @warn "Creating new figure" + visualization_callback.figure = GLMakie.Figure() + figure = visualization_callback.figure + if ndims == 2 + for v in 1:size(variable_names)[1] + GLMakie.heatmap(figure[makieLayoutHelper(v)...], plot_data.x, plot_data.y, plot_data.data[v]; + plot_arguments...) + end + else + for v in 1:size(variable_names)[1] + GLMakie.volume(figure[makieLayoutHelper(v)...], plot_data.data[v]; + algorithm = :mip, plot_arguments...) + end + end + GLMakie.display(figure) + else + if ndims == 2 + for v in 1:size(variable_names)[1] + GLMakie.heatmap!(figure[makieLayoutHelper(v)...], plot_data.x, plot_data.y, plot_data.data[v]; + plot_arguments...) + end + else + for v in 1:size(variable_names)[1] + GLMakie.volume!(figure[makieLayoutHelper(v)...], plot_data.data[v]; + algorithm = :mip, plot_arguments...) + end + end + end + + # TODO: show_mesh end - - GLMakie.display(visualization_callback.figure) - - # TODO: show_mesh -end -end # @muladd + end # @muladd + \ No newline at end of file From b1090ce229e138c50f89325c8073262a0a8cecc3 Mon Sep 17 00:00:00 2001 From: s6nistam Date: Thu, 12 Dec 2024 13:58:20 +0100 Subject: [PATCH 3/6] Visualization Callback with Makie working in 2D and 3D --- Project.toml | 2 + .../elixir_advection_amr_visualization.jl | 3 +- .../elixir_advection_amr_visualization.jl | 71 + src/callbacks_step/visualization.jl | 637 ++- src/visualization/types.jl | 135 +- src/visualization/utilities.jl | 3452 +++++++++-------- 6 files changed, 2461 insertions(+), 1839 deletions(-) create mode 100644 examples/tree_3d_dgsem/elixir_advection_amr_visualization.jl diff --git a/Project.toml b/Project.toml index 8a3a4772613..ba3f3117f47 100644 --- a/Project.toml +++ b/Project.toml @@ -15,6 +15,7 @@ Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -25,6 +26,7 @@ MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Preferences = "21216c6a-2e73-6563-6e65-726566657250" diff --git a/examples/tree_2d_dgsem/elixir_advection_amr_visualization.jl b/examples/tree_2d_dgsem/elixir_advection_amr_visualization.jl index 7b67b811177..07888f2f5b3 100644 --- a/examples/tree_2d_dgsem/elixir_advection_amr_visualization.jl +++ b/examples/tree_2d_dgsem/elixir_advection_amr_visualization.jl @@ -2,6 +2,7 @@ using OrdinaryDiffEq using Trixi using Plots +using GLMakie ############################################################################### # semidiscretization of the linear advection equation @@ -50,7 +51,7 @@ save_solution = SaveSolutionCallback(interval = 100, # Enable in-situ visualization with a new plot generated every 20 time steps # and additional plotting options passed as keyword arguments -visualization = VisualizationCallback(interval = 20, clims = (0, 1)) +visualization = VisualizationCallback(interval = 20, clims = (0, 1), plot_creator = Trixi.show_plot_makie) amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first), base_level = 3, diff --git a/examples/tree_3d_dgsem/elixir_advection_amr_visualization.jl b/examples/tree_3d_dgsem/elixir_advection_amr_visualization.jl new file mode 100644 index 00000000000..cb9b9c3cfca --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_advection_amr_visualization.jl @@ -0,0 +1,71 @@ + +using OrdinaryDiffEq +using Trixi +using Plots +using GLMakie + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +initial_condition = initial_condition_gauss +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + +coordinates_min = (-5.0, -5.0, -5.0) +coordinates_max = (5.0, 5.0, 5.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 30_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.3) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (entropy,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + + +visualization = VisualizationCallback(interval = 20, clims = (0, 1), plot_data_creator = Trixi.PlotData3D, plot_creator = Trixi.show_plot_makie) + +amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first), + base_level = 4, + med_level = 5, med_threshold = 0.1, + max_level = 6, max_threshold = 0.6) +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +stepsize_callback = StepsizeCallback(cfl = 1.2) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + visualization, + amr_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/callbacks_step/visualization.jl b/src/callbacks_step/visualization.jl index 3748f265eae..7d34e63f16e 100644 --- a/src/callbacks_step/visualization.jl +++ b/src/callbacks_step/visualization.jl @@ -3,339 +3,334 @@ # we need to opt-in explicitly. # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. @muladd begin - #! format: noindent - - mutable struct VisualizationCallback{SolutionVariables, VariableNames, PlotDataCreator, - PlotCreator} - interval::Int - solution_variables::SolutionVariables - variable_names::VariableNames - show_mesh::Bool - plot_data_creator::PlotDataCreator - plot_creator::PlotCreator - figure - plot_arguments::Dict{Symbol, Any} - end - - function Base.show(io::IO, - cb::DiscreteCallback{Condition, Affect!}) where {Condition, - Affect! <: - VisualizationCallback - } +#! format: noindent + +mutable struct VisualizationCallback{SolutionVariables, VariableNames, PlotDataCreator, + PlotCreator} + interval::Int + solution_variables::SolutionVariables + variable_names::VariableNames + show_mesh::Bool + plot_data_creator::PlotDataCreator + plot_creator::PlotCreator + figure + plot_arguments::Dict{Symbol, Any} +end + +function Base.show(io::IO, + cb::DiscreteCallback{Condition, Affect!}) where {Condition, + Affect! <: + VisualizationCallback + } + visualization_callback = cb.affect! + @unpack interval, plot_arguments, solution_variables, variable_names, show_mesh, plot_creator, plot_data_creator = visualization_callback + print(io, "VisualizationCallback(", + "interval=", interval, ", ", + "solution_variables=", solution_variables, ", ", + "variable_names=", variable_names, ", ", + "show_mesh=", show_mesh, ", ", + "plot_data_creator=", plot_data_creator, ", ", + "plot_creator=", plot_creator, ", ", + "plot_arguments=", plot_arguments, ")") +end + +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{Condition, Affect!}) where {Condition, + Affect! <: + VisualizationCallback + } + if get(io, :compact, false) + show(io, cb) + else visualization_callback = cb.affect! - @unpack interval, plot_arguments, solution_variables, variable_names, show_mesh, plot_creator, plot_data_creator = visualization_callback - print(io, "VisualizationCallback(", - "interval=", interval, ", ", - "solution_variables=", solution_variables, ", ", - "variable_names=", variable_names, ", ", - "show_mesh=", show_mesh, ", ", - "plot_data_creator=", plot_data_creator, ", ", - "plot_creator=", plot_creator, ", ", - "plot_arguments=", plot_arguments, ")") + + setup = [ + "interval" => visualization_callback.interval, + "plot arguments" => visualization_callback.plot_arguments, + "solution variables" => visualization_callback.solution_variables, + "variable names" => visualization_callback.variable_names, + "show mesh" => visualization_callback.show_mesh, + "plot creator" => visualization_callback.plot_creator, + "plot data creator" => visualization_callback.plot_data_creator + ] + summary_box(io, "VisualizationCallback", setup) end - - function Base.show(io::IO, ::MIME"text/plain", - cb::DiscreteCallback{Condition, Affect!}) where {Condition, - Affect! <: - VisualizationCallback - } - if get(io, :compact, false) - show(io, cb) - else - visualization_callback = cb.affect! - - setup = [ - "interval" => visualization_callback.interval, - "plot arguments" => visualization_callback.plot_arguments, - "solution variables" => visualization_callback.solution_variables, - "variable names" => visualization_callback.variable_names, - "show mesh" => visualization_callback.show_mesh, - "plot creator" => visualization_callback.plot_creator, - "plot data creator" => visualization_callback.plot_data_creator - ] - summary_box(io, "VisualizationCallback", setup) - end +end + +""" + VisualizationCallback(; interval=0, + solution_variables=cons2prim, + variable_names=[], + show_mesh=false, + plot_data_creator=PlotData2D, + plot_creator=show_plot, + plot_arguments...) + +Create a callback that visualizes results during a simulation, also known as *in-situ +visualization*. + +!!! warning "Experimental implementation" + This is an experimental feature and may change in any future releases. + +The `interval` specifies the number of time step iterations after which a new plot is generated. The +available variables to plot are configured with the `solution_variables` parameter, which acts the +same way as for the [`SaveSolutionCallback`](@ref). The variables to be actually plotted can be +selected by providing a single string or a list of strings to `variable_names`, and if `show_mesh` +is `true`, an additional plot with the mesh will be generated. + +To customize the generated figure, `plot_data_creator` allows to use different plot data types. With +`plot_creator` you can further specify an own function to visualize results, which must support the +same interface as the default implementation [`show_plot`](@ref). All remaining +keyword arguments are collected and passed as additional arguments to the plotting command. +""" +function VisualizationCallback(; interval = 0, + solution_variables = cons2prim, + variable_names = [], + show_mesh = false, + plot_data_creator = PlotData2D, + plot_creator = show_plot, + plot_arguments...) + mpi_isparallel() && error("this callback does not work in parallel yet") + + if variable_names isa String + variable_names = String[variable_names] end - - """ - VisualizationCallback(; interval=0, - solution_variables=cons2prim, - variable_names=[], - show_mesh=false, - plot_data_creator=PlotData2D, - plot_creator=show_plot, - plot_arguments...) - - Create a callback that visualizes results during a simulation, also known as *in-situ - visualization*. - - !!! warning "Experimental implementation" - This is an experimental feature and may change in any future releases. - - The `interval` specifies the number of time step iterations after which a new plot is generated. The - available variables to plot are configured with the `solution_variables` parameter, which acts the - same way as for the [`SaveSolutionCallback`](@ref). The variables to be actually plotted can be - selected by providing a single string or a list of strings to `variable_names`, and if `show_mesh` - is `true`, an additional plot with the mesh will be generated. - - To customize the generated figure, `plot_data_creator` allows to use different plot data types. With - `plot_creator` you can further specify an own function to visualize results, which must support the - same interface as the default implementation [`show_plot`](@ref). All remaining - keyword arguments are collected and passed as additional arguments to the plotting command. - """ - function VisualizationCallback(; interval = 0, - solution_variables = cons2prim, - variable_names = [], - show_mesh = false, - plot_data_creator = PlotData2D, - plot_creator = show_plot, - plot_arguments...) - mpi_isparallel() && error("this callback does not work in parallel yet") - - if variable_names isa String - variable_names = String[variable_names] - end - - visualization_callback = VisualizationCallback(interval, - solution_variables, variable_names, - show_mesh, - plot_data_creator, plot_creator, nothing, - Dict{Symbol, Any}(plot_arguments)) - - # Warn users if they create a visualization callback without having loaded the Plots package - # - # Note: This warning is added for convenience, as Plots is the only "officially" supported - # visualization package right now. However, in general nothing prevents anyone from using - # other packages such as Makie, Gadfly etc., given that appropriate `plot_creator`s are - # passed. This is also the reason why the visualization callback is not included via - # Requires.jl only when Plots is present. - # In the future, we should update/remove this warning if other plotting packages are - # starting to be used. - if !(:Plots in names(@__MODULE__, all = true) || :Makie in names(@__MODULE__, all = true)) - println(names(@__MODULE__, all = true)) - @warn "Package `Plots` or `Makie` required by `VisualizationCallback` to visualize results" - end - - DiscreteCallback(visualization_callback, visualization_callback, # the first one is the condition, the second the affect! - save_positions = (false, false), - initialize = initialize!) + + visualization_callback = VisualizationCallback(interval, + solution_variables, variable_names, + show_mesh, + plot_data_creator, plot_creator, nothing, + Dict{Symbol, Any}(plot_arguments)) + + # Warn users if they create a visualization callback without having loaded the Plots package + # + # Note: This warning is added for convenience, as Plots is the only "officially" supported + # visualization package right now. However, in general nothing prevents anyone from using + # other packages such as Makie, Gadfly etc., given that appropriate `plot_creator`s are + # passed. This is also the reason why the visualization callback is not included via + # Requires.jl only when Plots is present. + # In the future, we should update/remove this warning if other plotting packages are + # starting to be used. + if !(:Plots in names(@__MODULE__, all = true) || :Makie in names(@__MODULE__, all = true)) + println(names(@__MODULE__, all = true)) + @warn "Package `Plots` or `Makie` required by `VisualizationCallback` to visualize results" end - - function initialize!(cb::DiscreteCallback{Condition, Affect!}, u, t, - integrator) where {Condition, Affect! <: VisualizationCallback} - visualization_callback = cb.affect! - u_ode = integrator.u - semi = integrator.p - mesh, equations, solver, cache = mesh_equations_solver_cache(integrator.p) - if ndims(mesh) == 3 && visualization_callback.plot_data_creator == PlotData2D - visualization_callback.plot_data_creator = PlotData3D - end - visualization_callback(integrator) - return nothing + + DiscreteCallback(visualization_callback, visualization_callback, # the first one is the condition, the second the affect! + save_positions = (false, false), + initialize = initialize!) +end + +function initialize!(cb::DiscreteCallback{Condition, Affect!}, u, t, + integrator) where {Condition, Affect! <: VisualizationCallback} + visualization_callback = cb.affect! + u_ode = integrator.u + semi = integrator.p + mesh, equations, solver, cache = mesh_equations_solver_cache(integrator.p) + if ndims(mesh) == 3 && visualization_callback.plot_data_creator == PlotData2D + visualization_callback.plot_data_creator = PlotData3D end - - # this method is called to determine whether the callback should be activated - function (visualization_callback::VisualizationCallback)(u, t, integrator) - @unpack interval = visualization_callback - - # With error-based step size control, some steps can be rejected. Thus, - # `integrator.iter >= integrator.stats.naccept` - # (total #steps) (#accepted steps) - # We need to check the number of accepted steps since callbacks are not - # activated after a rejected step. - return interval > 0 && (integrator.stats.naccept % interval == 0 || - isfinished(integrator)) + visualization_callback(integrator) + return nothing +end + +# this method is called to determine whether the callback should be activated +function (visualization_callback::VisualizationCallback)(u, t, integrator) + @unpack interval = visualization_callback + + # With error-based step size control, some steps can be rejected. Thus, + # `integrator.iter >= integrator.stats.naccept` + # (total #steps) (#accepted steps) + # We need to check the number of accepted steps since callbacks are not + # activated after a rejected step. + return interval > 0 && (integrator.stats.naccept % interval == 0 || + isfinished(integrator)) +end + +# this method is called when the callback is activated +function (visualization_callback::VisualizationCallback)(integrator) + u_ode = integrator.u + semi = integrator.p + @unpack plot_arguments, solution_variables, variable_names, show_mesh, plot_data_creator, plot_creator, figure = visualization_callback + mesh, equations, solver, cache = mesh_equations_solver_cache(integrator.p) + n = Trixi.ndims(mesh) + + # Extract plot data + plot_data = plot_data_creator(u_ode, semi, solution_variables = solution_variables) + + # If variable names were not specified, plot everything + if isempty(variable_names) + variable_names = String[keys(plot_data)...] end - - # this method is called when the callback is activated - function (visualization_callback::VisualizationCallback)(integrator) - u_ode = integrator.u - semi = integrator.p - @unpack plot_arguments, solution_variables, variable_names, show_mesh, plot_data_creator, plot_creator, figure = visualization_callback - mesh, equations, solver, cache = mesh_equations_solver_cache(integrator.p) - n = Trixi.ndims(mesh) - - # Extract plot data - plot_data = plot_data_creator(u_ode, semi, solution_variables = solution_variables) - - # If variable names were not specified, plot everything - if isempty(variable_names) - variable_names = String[keys(plot_data)...] - end - - # Create plot - plot_creator(n, visualization_callback, plot_data, variable_names; - show_mesh = show_mesh, plot_arguments = plot_arguments, - time = integrator.t, timestep = integrator.stats.naccept, figure = figure) - - # avoid re-evaluating possible FSAL stages - u_modified!(integrator, false) - return nothing + + # Create plot + plot_creator(n, visualization_callback, plot_data, variable_names; + show_mesh = show_mesh, plot_arguments = plot_arguments, + time = integrator.t, timestep = integrator.stats.naccept, figure = figure) + + # avoid re-evaluating possible FSAL stages + u_modified!(integrator, false) + return nothing +end + +""" + show_plot(plot_data, variable_names; + show_mesh=true, plot_arguments=Dict{Symbol,Any}(), + time=nothing, timestep=nothing, figure=nothing) + +Visualize the plot data object provided in `plot_data` and display result, plotting only the +variables in `variable_names` and, optionally, the mesh (if `show_mesh` is `true`). Additionally, +`plot_arguments` will be unpacked and passed as keyword arguments to the `Plots.plot` command. + +This function is the default `plot_creator` argument for the [`VisualizationCallback`](@ref). +`time` and `timestep` are currently unused by this function. + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. + +See also: [`VisualizationCallback`](@ref), [`save_plot`](@ref) +""" +function show_plot(ndims, visualization_callback, plot_data, variable_names; + show_mesh = true, plot_arguments = Dict{Symbol, Any}(), + time = nothing, timestep = nothing, figure = nothing) + # Gather subplots + plots = [] + for v in variable_names + push!(plots, Plots.plot(plot_data[v]; plot_arguments...)) end - - """ - show_plot(plot_data, variable_names; - show_mesh=true, plot_arguments=Dict{Symbol,Any}(), - time=nothing, timestep=nothing, figure=nothing) - - Visualize the plot data object provided in `plot_data` and display result, plotting only the - variables in `variable_names` and, optionally, the mesh (if `show_mesh` is `true`). Additionally, - `plot_arguments` will be unpacked and passed as keyword arguments to the `Plots.plot` command. - - This function is the default `plot_creator` argument for the [`VisualizationCallback`](@ref). - `time` and `timestep` are currently unused by this function. - - !!! warning "Experimental implementation" - This is an experimental feature and may change in future releases. - - See also: [`VisualizationCallback`](@ref), [`save_plot`](@ref) - """ - function show_plot(ndims, visualization_callback, plot_data, variable_names; - show_mesh = true, plot_arguments = Dict{Symbol, Any}(), - time = nothing, timestep = nothing, figure = nothing) - # Gather subplots - plots = [] - for v in variable_names - push!(plots, Plots.plot(plot_data[v]; plot_arguments...)) - end - if show_mesh - push!(plots, Plots.plot(getmesh(plot_data); plot_arguments...)) - end - - # Note, for the visualization callback to work for general equation systems - # this layout construction would need to use the if-logic below. - # Currently, there is no use case for this so it is left here as a note. - # - # Determine layout - # if length(plots) <= 3 - # cols = length(plots) - # rows = 1 - # else - # cols = ceil(Int, sqrt(length(plots))) - # rows = div(length(plots), cols, RoundUp) - # end - # layout = (rows, cols) - - # Determine layout - cols = ceil(Int, sqrt(length(plots))) - rows = div(length(plots), cols, RoundUp) - layout = (rows, cols) - - # Show plot - display(Plots.plot(plots..., layout = layout)) + if show_mesh + push!(plots, Plots.plot(getmesh(plot_data); plot_arguments...)) end - - """ - save_plot(plot_data, variable_names; - show_mesh=true, plot_arguments=Dict{Symbol,Any}(), - time=nothing, timestep=nothing) - - Visualize the plot data object provided in `plot_data` and save result as a PNG file in the `out` - directory, plotting only the variables in `variable_names` and, optionally, the mesh (if `show_mesh` - is `true`). Additionally, `plot_arguments` will be unpacked and passed as keyword arguments to the - `Plots.plot` command. - - The `timestep` is used in the filename. `time` is currently unused by this function. - - !!! warning "Experimental implementation" - This is an experimental feature and may change in future releases. - - See also: [`VisualizationCallback`](@ref), [`show_plot`](@ref) - """ - function save_plot(plot_data, variable_names; - show_mesh = true, plot_arguments = Dict{Symbol, Any}(), - time = nothing, timestep = nothing) - # Gather subplots - plots = [] - for v in variable_names - push!(plots, Plots.plot(plot_data[v]; plot_arguments...)) - end - if show_mesh - push!(plots, Plots.plot(getmesh(plot_data); plot_arguments...)) - end - - # Determine layout - cols = ceil(Int, sqrt(length(plots))) - rows = div(length(plots), cols, RoundUp) - layout = (rows, cols) - - # Create plot - Plots.plot(plots..., layout = layout) - - # Determine filename and save plot - filename = joinpath("out", @sprintf("solution_%09d.png", timestep)) - Plots.savefig(filename) + + # Note, for the visualization callback to work for general equation systems + # this layout construction would need to use the if-logic below. + # Currently, there is no use case for this so it is left here as a note. + # + # Determine layout + # if length(plots) <= 3 + # cols = length(plots) + # rows = 1 + # else + # cols = ceil(Int, sqrt(length(plots))) + # rows = div(length(plots), cols, RoundUp) + # end + # layout = (rows, cols) + + # Determine layout + cols = ceil(Int, sqrt(length(plots))) + rows = div(length(plots), cols, RoundUp) + layout = (rows, cols) + + # Show plot + display(Plots.plot(plots..., layout = layout)) +end + +""" + save_plot(plot_data, variable_names; + show_mesh=true, plot_arguments=Dict{Symbol,Any}(), + time=nothing, timestep=nothing) + +Visualize the plot data object provided in `plot_data` and save result as a PNG file in the `out` +directory, plotting only the variables in `variable_names` and, optionally, the mesh (if `show_mesh` +is `true`). Additionally, `plot_arguments` will be unpacked and passed as keyword arguments to the +`Plots.plot` command. + +The `timestep` is used in the filename. `time` is currently unused by this function. + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. + +See also: [`VisualizationCallback`](@ref), [`show_plot`](@ref) +""" +function save_plot(plot_data, variable_names; + show_mesh = true, plot_arguments = Dict{Symbol, Any}(), + time = nothing, timestep = nothing) + # Gather subplots + plots = [] + for v in variable_names + push!(plots, Plots.plot(plot_data[v]; plot_arguments...)) end - - #converts a single int into a tuple of ints, to get a square arrangement for example f(1) = (1,1) f(2) = (2,1) f(3) = (2,2) f(4) = (1,2) - function makieLayoutHelper(n) - if n == 1 - return (1,1) - end - t = intTo2DInt(n-1) - if t[1] == 1 - return (t[2] + 1, 1) - elseif t[1] > t[2] - return (t[1], t[2] + 1) - elseif t[2] >= t[1] - return (t[1] - 1, t[2]) - end + if show_mesh + push!(plots, Plots.plot(getmesh(plot_data); plot_arguments...)) + end + + # Determine layout + cols = ceil(Int, sqrt(length(plots))) + rows = div(length(plots), cols, RoundUp) + layout = (rows, cols) + + # Create plot + Plots.plot(plots..., layout = layout) + + # Determine filename and save plot + filename = joinpath("out", @sprintf("solution_%09d.png", timestep)) + Plots.savefig(filename) +end + +#converts a single int into a tuple of ints, to get a square arrangement for example f(1) = (1,1) f(2) = (2,1) f(3) = (2,2) f(4) = (1,2) +function makieLayoutHelper(n) + if n == 1 + return (1,1) end - - - """ - show_plot_makie(plot_data, variable_names; - show_mesh=true, plot_arguments=Dict{Symbol,Any}(), - time=nothing, timestep=nothing, figure=nothing) - - Visualize the plot data object provided in `plot_data` and display result using Makie. - Only the variables in `variable_names` are plotted and, optionally, the mesh (if - `show_mesh` is `true`). Additionally, `plot_arguments` will be unpacked and passed as - keyword arguments to the `Plots.plot` command. - - This function is the default `plot_creator` argument for the [`VisualizationCallback`](@ref). - `time` and `timestep` are currently unused by this function. - - !!! warning "Experimental implementation" - This is an experimental feature and may change in future releases. - - See also: [`VisualizationCallback`](@ref), [`save_plot`](@ref) - """ - function show_plot_makie(ndims, visualization_callback, plot_data, variable_names; - show_mesh = true, plot_arguments = Dict{Symbol, Any}(), - time = nothing, timestep = nothing, figure = nothing) - - if figure === nothing - @warn "Creating new figure" - visualization_callback.figure = GLMakie.Figure() - figure = visualization_callback.figure - if ndims == 2 - for v in 1:size(variable_names)[1] - GLMakie.heatmap(figure[makieLayoutHelper(v)...], plot_data.x, plot_data.y, plot_data.data[v]; - plot_arguments...) - end - else - for v in 1:size(variable_names)[1] - GLMakie.volume(figure[makieLayoutHelper(v)...], plot_data.data[v]; - algorithm = :mip, plot_arguments...) - end + t = makieLayoutHelper(n-1) + if t[1] == 1 + return (t[2] + 1, 1) + elseif t[1] > t[2] + return (t[1], t[2] + 1) + elseif t[2] >= t[1] + return (t[1] - 1, t[2]) + end +end + + +""" + show_plot_makie(plot_data, variable_names; + show_mesh=true, plot_arguments=Dict{Symbol,Any}(), + time=nothing, timestep=nothing, figure=nothing) + +Visualize the plot data object provided in `plot_data` and display result using Makie. +Only the variables in `variable_names` are plotted and, optionally, the mesh (if +`show_mesh` is `true`). Additionally, `plot_arguments` will be unpacked and passed as +keyword arguments to the `Plots.plot` command. + +This function is the default `plot_creator` argument for the [`VisualizationCallback`](@ref). +`time` and `timestep` are currently unused by this function. + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. + +See also: [`VisualizationCallback`](@ref), [`save_plot`](@ref) +""" +function show_plot_makie(ndims, visualization_callback, plot_data, variable_names; + show_mesh = true, plot_arguments = Dict{Symbol, Any}(), + time = nothing, timestep = nothing, figure = nothing) + + if figure === nothing + @warn "Creating new figure" + visualization_callback.figure = GLMakie.Figure() + figure = visualization_callback.figure + if ndims == 2 + for v in 1:size(variable_names)[1] + GLMakie.heatmap(figure[makieLayoutHelper(v)...], plot_data.x, plot_data.y, plot_data.data[v]) + end + else + for v in 1:size(variable_names)[1] + GLMakie.volume(figure[makieLayoutHelper(v)...], plot_data.data[v]) + end + end + GLMakie.display(figure) + else + if ndims == 2 + for v in 1:size(variable_names)[1] + GLMakie.heatmap!(figure[makieLayoutHelper(v)...], plot_data.x, plot_data.y, plot_data.data[v]) end - GLMakie.display(figure) else - if ndims == 2 - for v in 1:size(variable_names)[1] - GLMakie.heatmap!(figure[makieLayoutHelper(v)...], plot_data.x, plot_data.y, plot_data.data[v]; - plot_arguments...) - end - else - for v in 1:size(variable_names)[1] - GLMakie.volume!(figure[makieLayoutHelper(v)...], plot_data.data[v]; - algorithm = :mip, plot_arguments...) - end + for v in 1:size(variable_names)[1] + GLMakie.volume!(figure[makieLayoutHelper(v)...], plot_data.data[v]) end end - - # TODO: show_mesh end - end # @muladd - \ No newline at end of file + + # TODO: show_mesh +end +end # @muladd diff --git a/src/visualization/types.jl b/src/visualization/types.jl index b294ce25607..bf3f4c0ad59 100644 --- a/src/visualization/types.jl +++ b/src/visualization/types.jl @@ -58,6 +58,42 @@ end Base.eltype(pd::AbstractPlotData) = Pair{String, PlotDataSeries{typeof(pd)}} +""" + PlotData3D + +Holds all relevant data for creating 3D plots of multiple solution variables and to visualize the +mesh. + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. +""" +struct PlotData3DCartesian{Coordinates, Data, VariableNames, Vertices} <: + AbstractPlotData{3} + x::Coordinates + y::Coordinates + z::Coordinates + data::Data + variable_names::VariableNames + mesh_vertices_x::Vertices + mesh_vertices_y::Vertices + mesh_vertices_z::Vertices + orientation_x::Int + orientation_y::Int + orientation_z::Int +end + +# Show only a truncated output for convenience (the full data does not make sense) +function Base.show(io::IO, pd::PlotData3DCartesian) + @nospecialize pd # reduce precompilation time + + print(io, "PlotData3DCartesian{", + typeof(pd.x), ",", + typeof(pd.data), ",", + typeof(pd.variable_names), ",", + typeof(pd.mesh_vertices_x), + "}(, , , , , , , )") +end + """ PlotData2D @@ -183,6 +219,103 @@ Extract grid lines from `pd` for plotting with `Plots.plot`. """ getmesh(pd::AbstractPlotData) = PlotMesh(pd) +""" + PlotData3D(u, semi [or mesh, equations, solver, cache]; + solution_variables=nothing, + grid_lines=true, max_supported_level=11, nvisnodes=nothing, + slice=:xy, point=(0.0, 0.0, 0.0)) + +Create a new `PlotData3D` object that can be used for visualizing 2D/3D DGSEM solution data array +`u` with `Plots.jl`. All relevant geometrical information is extracted from the semidiscretization +`semi`. By default, the primitive variables (if existent) or the conservative variables (otherwise) +from the solution are used for plotting. This can be changed by passing an appropriate conversion +function to `solution_variables`. + +If `grid_lines` is `true`, also extract grid vertices for visualizing the mesh. The output +resolution is indirectly set via `max_supported_level`: all data is interpolated to +`2^max_supported_level` uniformly distributed points in each spatial direction, also setting the +maximum allowed refinement level in the solution. `nvisnodes` specifies the number of visualization +nodes to be used. If it is `nothing`, twice the number of solution DG nodes are used for +visualization, and if set to `0`, exactly the number of nodes in the DG elements are used. + +When visualizing data from a three-dimensional simulation, a 2D slice is extracted for plotting. +`slice` specifies the plane that is being sliced and may be `:xy`, `:xz`, or `:yz`. +The slice position is specified by a `point` that lies on it, which defaults to `(0.0, 0.0, 0.0)`. +Both of these values are ignored when visualizing 2D data. + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. + +# Examples +```julia +julia> using Trixi, Plots + +julia> trixi_include(default_example()) +[...] + +julia> pd = PlotData3D(sol) +PlotData3D(...) + +julia> plot(pd) # To plot all available variables + +julia> plot(pd["scalar"]) # To plot only a single variable + +julia> plot!(getmesh(pd)) # To add grid lines to the plot +``` +""" + +function PlotData3D(sol::TrixiODESolution; kwargs...) + PlotData3D(sol.u[end], sol.prob.p; kwargs...) +end + +function PlotData3D(u_ode, semi; kwargs...) + PlotData3D(wrap_array_native(u_ode, semi), + mesh_equations_solver_cache(semi)...; + kwargs...) +end + +# Create a PlotData3DCartesian object for TreeMeshes on default. +function PlotData3D(u, mesh::TreeMesh, equations, solver, cache; kwargs...) + PlotData3DCartesian(u, mesh::TreeMesh, equations, solver, cache; kwargs...) +end + +# Create a PlotData3DCartesian for a TreeMesh. +function PlotData3DCartesian(u, mesh::TreeMesh, equations, solver, cache; + solution_variables = nothing, + grid_lines = true, max_supported_level = 11, + nvisnodes = nothing) + @assert ndims(mesh)==3 "unsupported number of dimensions $ndims (must be 3)" + solution_variables_ = digest_solution_variables(equations, solution_variables) + + # Extract mesh info + center_level_0 = mesh.tree.center_level_0 + length_level_0 = mesh.tree.length_level_0 + leaf_cell_ids = leaf_cells(mesh.tree) + coordinates = mesh.tree.coordinates[:, leaf_cell_ids] + levels = mesh.tree.levels[leaf_cell_ids] + + unstructured_data = get_unstructured_data(u, solution_variables_, mesh, equations, + solver, cache) + x, y, z, data, mesh_vertices_x, mesh_vertices_y, mesh_vertices_z = get_data_3d(center_level_0, + length_level_0, + coordinates, + levels, + unstructured_data, + nnodes(solver), + grid_lines, + max_supported_level, + nvisnodes) + variable_names = SVector(varnames(solution_variables_, equations)) + + orientation_x = 1 + orientation_y = 2 + orientation_z = 3 + + return PlotData3DCartesian(x, y, z, data, variable_names, mesh_vertices_x, + mesh_vertices_y, mesh_vertices_z, + orientation_x, orientation_y, orientation_z) +end + """ PlotData2D(u, semi [or mesh, equations, solver, cache]; solution_variables=nothing, @@ -719,4 +852,4 @@ function PlotData1D(cb::DiscreteCallback{<:Any, <:TimeSeriesCallback}, point_id::Integer) return PlotData1D(cb.affect!, point_id) end -end # @muladd +end # @muladd \ No newline at end of file diff --git a/src/visualization/utilities.jl b/src/visualization/utilities.jl index 1f843c6a9d2..324219259b4 100644 --- a/src/visualization/utilities.jl +++ b/src/visualization/utilities.jl @@ -3,1576 +3,1996 @@ # we need to opt-in explicitly. # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. @muladd begin -#! format: noindent - -@inline num_faces(elem::Tri) = 3 -@inline num_faces(elem::Quad) = 4 - -# compute_triangle_area(tri) -# -# Computes the area of a triangle given `tri`, which is a tuple of three points (vectors), -# using the [Shoelace_formula](https://en.wikipedia.org/wiki/Shoelace_formula). -function compute_triangle_area(tri) - A, B, C = tri - return 0.5f0 * (A[1] * (B[2] - C[2]) + B[1] * (C[2] - A[2]) + C[1] * (A[2] - B[2])) -end - -# reference_plotting_triangulation(reference_plotting_coordinates) -# -# Computes a triangulation of the points in `reference_plotting_coordinates`, which is a tuple containing -# vectors of plotting points on the reference element (e.g., reference_plotting_coordinates = (r,s)). -# The reference element is assumed to be [-1,1]^d. -# -# This function returns `t` which is a `3 x N_tri` Matrix{Int} containing indices of triangles in the -# triangulation of the plotting points, with zero-volume triangles removed. -# -# For example, r[t[1, i]] returns the first reference coordinate of the 1st point on the ith triangle. -function reference_plotting_triangulation(reference_plotting_coordinates, - tol = 50 * eps()) - # on-the-fly triangulation of plotting nodes on the reference element - tri_in = Triangulate.TriangulateIO() - tri_in.pointlist = permutedims(hcat(reference_plotting_coordinates...)) - tri_out, _ = Triangulate.triangulate("Q", tri_in) - triangles = tri_out.trianglelist - - # filter out sliver triangles - has_volume = fill(true, size(triangles, 2)) - for i in axes(triangles, 2) - ids = @view triangles[:, i] - x_points = @view tri_out.pointlist[1, ids] - y_points = @view tri_out.pointlist[2, ids] - area = compute_triangle_area(zip(x_points, y_points)) - if abs(area) < tol - has_volume[i] = false + #! format: noindent + + @inline num_faces(elem::Tri) = 3 + @inline num_faces(elem::Quad) = 4 + + # compute_triangle_area(tri) + # + # Computes the area of a triangle given `tri`, which is a tuple of three points (vectors), + # using the [Shoelace_formula](https://en.wikipedia.org/wiki/Shoelace_formula). + function compute_triangle_area(tri) + A, B, C = tri + return 0.5f0 * (A[1] * (B[2] - C[2]) + B[1] * (C[2] - A[2]) + C[1] * (A[2] - B[2])) + end + + # reference_plotting_triangulation(reference_plotting_coordinates) + # + # Computes a triangulation of the points in `reference_plotting_coordinates`, which is a tuple containing + # vectors of plotting points on the reference element (e.g., reference_plotting_coordinates = (r,s)). + # The reference element is assumed to be [-1,1]^d. + # + # This function returns `t` which is a `3 x N_tri` Matrix{Int} containing indices of triangles in the + # triangulation of the plotting points, with zero-volume triangles removed. + # + # For example, r[t[1, i]] returns the first reference coordinate of the 1st point on the ith triangle. + function reference_plotting_triangulation(reference_plotting_coordinates, + tol = 50 * eps()) + # on-the-fly triangulation of plotting nodes on the reference element + tri_in = Triangulate.TriangulateIO() + tri_in.pointlist = permutedims(hcat(reference_plotting_coordinates...)) + tri_out, _ = Triangulate.triangulate("Q", tri_in) + triangles = tri_out.trianglelist + + # filter out sliver triangles + has_volume = fill(true, size(triangles, 2)) + for i in axes(triangles, 2) + ids = @view triangles[:, i] + x_points = @view tri_out.pointlist[1, ids] + y_points = @view tri_out.pointlist[2, ids] + area = compute_triangle_area(zip(x_points, y_points)) + if abs(area) < tol + has_volume[i] = false + end end + return permutedims(triangles[:, findall(has_volume)]) end - return permutedims(triangles[:, findall(has_volume)]) -end - -# This function is used to avoid type instabilities when calling `digest_solution_variables`. -function transform_to_solution_variables!(u, solution_variables, equations) - for (i, u_i) in enumerate(u) - u[i] = solution_variables(u_i, equations) - end -end - -# global_plotting_triangulation_triplot(u_plot, rst_plot, xyz_plot) -# -# Returns (plotting_coordinates_x, plotting_coordinates_y, ..., plotting_values, plotting_triangulation). -# Output can be used with TriplotRecipes.DGTriPseudocolor(...). -# -# Inputs: -# - xyz_plot = plotting points (tuple of matrices of size (Nplot, K)) -# - u_plot = matrix of size (Nplot, K) representing solution to plot. -# - t = triangulation of reference plotting points -function global_plotting_triangulation_triplot(xyz_plot, u_plot, t) - @assert size(first(xyz_plot), 1)==size(u_plot, 1) "Row dimension of u_plot does not match row dimension of xyz_plot" - - # build discontinuous data on plotting triangular mesh - num_plotting_points, num_elements = size(u_plot) - num_reference_plotting_triangles = size(t, 1) - num_plotting_elements_total = num_reference_plotting_triangles * num_elements - - # each column of `tp` corresponds to a vertex of a plotting triangle - tp = zeros(Int32, 3, num_plotting_elements_total) - zp = similar(tp, eltype(u_plot)) - for e in 1:num_elements - for i in 1:num_reference_plotting_triangles - tp[:, i + (e - 1) * num_reference_plotting_triangles] .= @views t[i, :] .+ - (e - 1) * - num_plotting_points - zp[:, i + (e - 1) * num_reference_plotting_triangles] .= @views u_plot[t[i, - :], - e] + + # This function is used to avoid type instabilities when calling `digest_solution_variables`. + function transform_to_solution_variables!(u, solution_variables, equations) + for (i, u_i) in enumerate(u) + u[i] = solution_variables(u_i, equations) end end - return vec.(xyz_plot)..., zp, tp -end - -function get_face_node_indices(r, s, dg::DGSEM, tol = 100 * eps()) - face_1 = findall(@. abs(s + 1) < tol) - face_2 = findall(@. abs(r - 1) < tol) - face_3 = findall(@. abs(s - 1) < tol) - face_4 = findall(@. abs(r + 1) < tol) - Fmask = hcat(face_1, face_2, face_3, face_4) - return Fmask -end - -# dispatch on semi -function mesh_plotting_wireframe(u, semi) - mesh_plotting_wireframe(u, mesh_equations_solver_cache(semi)...) -end - -# mesh_plotting_wireframe(u, mesh, equations, dg::DGMulti, cache; num_plotting_pts=25) -# -# Generates data for plotting a mesh wireframe given StartUpDG data types. -# Returns (plotting_coordinates_x, plotting_coordinates_y, nothing) for a 2D mesh wireframe. -function mesh_plotting_wireframe(u::StructArray, mesh, equations, dg::DGMulti, cache; - nvisnodes = 2 * nnodes(dg)) - @unpack md = mesh - rd = dg.basis - - # Construct 1D plotting interpolation matrix `Vp1D` for a single face - @unpack N, Fmask = rd - num_face_points = length(Fmask) ÷ num_faces(rd.element_type) - vandermonde_matrix_1D = StartUpDG.vandermonde(Line(), N, - StartUpDG.nodes(Line(), - num_face_points - 1)) - rplot = LinRange(-1, 1, nvisnodes) - Vp1D = StartUpDG.vandermonde(Line(), N, rplot) / vandermonde_matrix_1D - - num_faces_total = num_faces(rd.element_type) * md.num_elements - xf, yf = map(x -> reshape(view(x, Fmask, :), num_face_points, num_faces_total), - md.xyz) - uf = similar(u, size(xf)) - apply_to_each_field((out, x) -> out .= reshape(view(x, Fmask, :), num_face_points, - num_faces_total), uf, u) - - num_face_plotting_points = size(Vp1D, 1) - x_mesh, y_mesh = ntuple(_ -> zeros(num_face_plotting_points, num_faces_total), 2) - u_mesh = similar(u, (num_face_plotting_points, num_faces_total)) - for f in 1:num_faces_total - mul!(view(x_mesh, :, f), Vp1D, view(xf, :, f)) - mul!(view(y_mesh, :, f), Vp1D, view(yf, :, f)) - apply_to_each_field(mul_by!(Vp1D), view(u_mesh, :, f), view(uf, :, f)) - end - - return x_mesh, y_mesh, u_mesh -end - -function mesh_plotting_wireframe(u::StructArray, mesh, equations, dg::DGSEM, cache; - nvisnodes = 2 * nnodes(dg)) - - # build nodes on reference element (seems to be the right ordering) - r, s = reference_node_coordinates_2d(dg) - - # extract node coordinates - uEltype = eltype(first(u)) - nvars = nvariables(equations) - n_nodes_2d = nnodes(dg)^ndims(mesh) - n_elements = nelements(dg, cache) - x = reshape(view(cache.elements.node_coordinates, 1, :, :, :), n_nodes_2d, - n_elements) - y = reshape(view(cache.elements.node_coordinates, 2, :, :, :), n_nodes_2d, - n_elements) - - # extract indices of local face nodes for wireframe plotting - Fmask = get_face_node_indices(r, s, dg) - plotting_interp_matrix1D = face_plotting_interpolation_matrix(dg; - nvisnodes = nvisnodes) - - # These 5 lines extract the face values on each element from the arrays x,y,sol_to_plot. - # The resulting arrays are then reshaped so that xf, yf, sol_f are Matrix types of size - # (Number of face plotting nodes) x (Number of faces). - function face_first_reshape(x, num_nodes_1D, num_nodes, num_elements) - num_reference_faces = 2 * ndims(mesh) - xf = view(reshape(x, num_nodes, num_elements), vec(Fmask), :) - return reshape(xf, num_nodes_1D, num_elements * num_reference_faces) - end - function reshape_and_interpolate(x) - plotting_interp_matrix1D * - face_first_reshape(x, nnodes(dg), n_nodes_2d, n_elements) - end - xfp, yfp = map(reshape_and_interpolate, (x, y)) - ufp = StructArray{SVector{nvars, uEltype}}(map(reshape_and_interpolate, - StructArrays.components(u))) - - return xfp, yfp, ufp -end - -function mesh_plotting_wireframe(u::ScalarData, mesh, equations, dg::DGSEM, cache; - nvisnodes = 2 * nnodes(dg)) - - # build nodes on reference element (seems to be the right ordering) - r, s = reference_node_coordinates_2d(dg) - - # extract node coordinates - n_nodes_2d = nnodes(dg)^ndims(mesh) - n_elements = nelements(dg, cache) - x = reshape(view(cache.elements.node_coordinates, 1, :, :, :), n_nodes_2d, - n_elements) - y = reshape(view(cache.elements.node_coordinates, 2, :, :, :), n_nodes_2d, - n_elements) - - # extract indices of local face nodes for wireframe plotting - Fmask = get_face_node_indices(r, s, dg) - plotting_interp_matrix1D = face_plotting_interpolation_matrix(dg; - nvisnodes = nvisnodes) - - # These 5 lines extract the face values on each element from the arrays x,y,sol_to_plot. - # The resulting arrays are then reshaped so that xf, yf, sol_f are Matrix types of size - # (Number of face plotting nodes) x (Number of faces). - function face_first_reshape(x, num_nodes_1D, num_nodes, num_elements) - num_reference_faces = 2 * ndims(mesh) - xf = view(reshape(x, num_nodes, num_elements), vec(Fmask), :) - return reshape(xf, num_nodes_1D, num_elements * num_reference_faces) - end - function reshape_and_interpolate(x) - plotting_interp_matrix1D * - face_first_reshape(x, nnodes(dg), n_nodes_2d, n_elements) - end - xfp, yfp, ufp = map(reshape_and_interpolate, (x, y, u.data)) - - return xfp, yfp, ufp -end - -function mesh_plotting_wireframe(u::ScalarData, mesh, equations, dg::DGMulti, cache; - nvisnodes = 2 * nnodes(dg)) - @unpack md = mesh - rd = dg.basis - - # Construct 1D plotting interpolation matrix `Vp1D` for a single face - @unpack N, Fmask = rd - vandermonde_matrix_1D = StartUpDG.vandermonde(Line(), N, StartUpDG.nodes(Line(), N)) - rplot = LinRange(-1, 1, nvisnodes) - Vp1D = StartUpDG.vandermonde(Line(), N, rplot) / vandermonde_matrix_1D - - num_face_points = N + 1 - num_faces_total = num_faces(rd.element_type) * md.num_elements - xf, yf, uf = map(x -> reshape(view(x, Fmask, :), num_face_points, num_faces_total), - (md.xyz..., u.data)) - - num_face_plotting_points = size(Vp1D, 1) - x_mesh, y_mesh = ntuple(_ -> zeros(num_face_plotting_points, num_faces_total), 2) - u_mesh = similar(u.data, (num_face_plotting_points, num_faces_total)) - for f in 1:num_faces_total - mul!(view(x_mesh, :, f), Vp1D, view(xf, :, f)) - mul!(view(y_mesh, :, f), Vp1D, view(yf, :, f)) - mul!(view(u_mesh, :, f), Vp1D, view(uf, :, f)) - end - - return x_mesh, y_mesh, u_mesh -end - -# These methods are used internally to set the default value of the solution variables: -# - If a `cons2prim` for the given `equations` exists, use it -# - Otherwise, use `cons2cons`, which is defined for all systems of equations -digest_solution_variables(equations, solution_variables) = solution_variables -function digest_solution_variables(equations, solution_variables::Nothing) - if hasmethod(cons2prim, Tuple{AbstractVector, typeof(equations)}) - return cons2prim - else - return cons2cons + + # global_plotting_triangulation_triplot(u_plot, rst_plot, xyz_plot) + # + # Returns (plotting_coordinates_x, plotting_coordinates_y, ..., plotting_values, plotting_triangulation). + # Output can be used with TriplotRecipes.DGTriPseudocolor(...). + # + # Inputs: + # - xyz_plot = plotting points (tuple of matrices of size (Nplot, K)) + # - u_plot = matrix of size (Nplot, K) representing solution to plot. + # - t = triangulation of reference plotting points + function global_plotting_triangulation_triplot(xyz_plot, u_plot, t) + @assert size(first(xyz_plot), 1)==size(u_plot, 1) "Row dimension of u_plot does not match row dimension of xyz_plot" + + # build discontinuous data on plotting triangular mesh + num_plotting_points, num_elements = size(u_plot) + num_reference_plotting_triangles = size(t, 1) + num_plotting_elements_total = num_reference_plotting_triangles * num_elements + + # each column of `tp` corresponds to a vertex of a plotting triangle + tp = zeros(Int32, 3, num_plotting_elements_total) + zp = similar(tp, eltype(u_plot)) + for e in 1:num_elements + for i in 1:num_reference_plotting_triangles + tp[:, i + (e - 1) * num_reference_plotting_triangles] .= @views t[i, :] .+ + (e - 1) * + num_plotting_points + zp[:, i + (e - 1) * num_reference_plotting_triangles] .= @views u_plot[t[i, + :], + e] + end + end + return vec.(xyz_plot)..., zp, tp + end + + function get_face_node_indices(r, s, dg::DGSEM, tol = 100 * eps()) + face_1 = findall(@. abs(s + 1) < tol) + face_2 = findall(@. abs(r - 1) < tol) + face_3 = findall(@. abs(s - 1) < tol) + face_4 = findall(@. abs(r + 1) < tol) + Fmask = hcat(face_1, face_2, face_3, face_4) + return Fmask + end + + # dispatch on semi + function mesh_plotting_wireframe(u, semi) + mesh_plotting_wireframe(u, mesh_equations_solver_cache(semi)...) + end + + # mesh_plotting_wireframe(u, mesh, equations, dg::DGMulti, cache; num_plotting_pts=25) + # + # Generates data for plotting a mesh wireframe given StartUpDG data types. + # Returns (plotting_coordinates_x, plotting_coordinates_y, nothing) for a 2D mesh wireframe. + function mesh_plotting_wireframe(u::StructArray, mesh, equations, dg::DGMulti, cache; + nvisnodes = 2 * nnodes(dg)) + @unpack md = mesh + rd = dg.basis + + # Construct 1D plotting interpolation matrix `Vp1D` for a single face + @unpack N, Fmask = rd + num_face_points = length(Fmask) ÷ num_faces(rd.element_type) + vandermonde_matrix_1D = StartUpDG.vandermonde(Line(), N, + StartUpDG.nodes(Line(), + num_face_points - 1)) + rplot = LinRange(-1, 1, nvisnodes) + Vp1D = StartUpDG.vandermonde(Line(), N, rplot) / vandermonde_matrix_1D + + num_faces_total = num_faces(rd.element_type) * md.num_elements + xf, yf = map(x -> reshape(view(x, Fmask, :), num_face_points, num_faces_total), + md.xyz) + uf = similar(u, size(xf)) + apply_to_each_field((out, x) -> out .= reshape(view(x, Fmask, :), num_face_points, + num_faces_total), uf, u) + + num_face_plotting_points = size(Vp1D, 1) + x_mesh, y_mesh = ntuple(_ -> zeros(num_face_plotting_points, num_faces_total), 2) + u_mesh = similar(u, (num_face_plotting_points, num_faces_total)) + for f in 1:num_faces_total + mul!(view(x_mesh, :, f), Vp1D, view(xf, :, f)) + mul!(view(y_mesh, :, f), Vp1D, view(yf, :, f)) + apply_to_each_field(mul_by!(Vp1D), view(u_mesh, :, f), view(uf, :, f)) + end + + return x_mesh, y_mesh, u_mesh + end + + function mesh_plotting_wireframe(u::StructArray, mesh, equations, dg::DGSEM, cache; + nvisnodes = 2 * nnodes(dg)) + + # build nodes on reference element (seems to be the right ordering) + r, s = reference_node_coordinates_2d(dg) + + # extract node coordinates + uEltype = eltype(first(u)) + nvars = nvariables(equations) + n_nodes_2d = nnodes(dg)^ndims(mesh) + n_elements = nelements(dg, cache) + x = reshape(view(cache.elements.node_coordinates, 1, :, :, :), n_nodes_2d, + n_elements) + y = reshape(view(cache.elements.node_coordinates, 2, :, :, :), n_nodes_2d, + n_elements) + + # extract indices of local face nodes for wireframe plotting + Fmask = get_face_node_indices(r, s, dg) + plotting_interp_matrix1D = face_plotting_interpolation_matrix(dg; + nvisnodes = nvisnodes) + + # These 5 lines extract the face values on each element from the arrays x,y,sol_to_plot. + # The resulting arrays are then reshaped so that xf, yf, sol_f are Matrix types of size + # (Number of face plotting nodes) x (Number of faces). + function face_first_reshape(x, num_nodes_1D, num_nodes, num_elements) + num_reference_faces = 2 * ndims(mesh) + xf = view(reshape(x, num_nodes, num_elements), vec(Fmask), :) + return reshape(xf, num_nodes_1D, num_elements * num_reference_faces) + end + function reshape_and_interpolate(x) + plotting_interp_matrix1D * + face_first_reshape(x, nnodes(dg), n_nodes_2d, n_elements) + end + xfp, yfp = map(reshape_and_interpolate, (x, y)) + ufp = StructArray{SVector{nvars, uEltype}}(map(reshape_and_interpolate, + StructArrays.components(u))) + + return xfp, yfp, ufp + end + + function mesh_plotting_wireframe(u::ScalarData, mesh, equations, dg::DGSEM, cache; + nvisnodes = 2 * nnodes(dg)) + + # build nodes on reference element (seems to be the right ordering) + r, s = reference_node_coordinates_2d(dg) + + # extract node coordinates + n_nodes_2d = nnodes(dg)^ndims(mesh) + n_elements = nelements(dg, cache) + x = reshape(view(cache.elements.node_coordinates, 1, :, :, :), n_nodes_2d, + n_elements) + y = reshape(view(cache.elements.node_coordinates, 2, :, :, :), n_nodes_2d, + n_elements) + + # extract indices of local face nodes for wireframe plotting + Fmask = get_face_node_indices(r, s, dg) + plotting_interp_matrix1D = face_plotting_interpolation_matrix(dg; + nvisnodes = nvisnodes) + + # These 5 lines extract the face values on each element from the arrays x,y,sol_to_plot. + # The resulting arrays are then reshaped so that xf, yf, sol_f are Matrix types of size + # (Number of face plotting nodes) x (Number of faces). + function face_first_reshape(x, num_nodes_1D, num_nodes, num_elements) + num_reference_faces = 2 * ndims(mesh) + xf = view(reshape(x, num_nodes, num_elements), vec(Fmask), :) + return reshape(xf, num_nodes_1D, num_elements * num_reference_faces) + end + function reshape_and_interpolate(x) + plotting_interp_matrix1D * + face_first_reshape(x, nnodes(dg), n_nodes_2d, n_elements) + end + xfp, yfp, ufp = map(reshape_and_interpolate, (x, y, u.data)) + + return xfp, yfp, ufp + end + + function mesh_plotting_wireframe(u::ScalarData, mesh, equations, dg::DGMulti, cache; + nvisnodes = 2 * nnodes(dg)) + @unpack md = mesh + rd = dg.basis + + # Construct 1D plotting interpolation matrix `Vp1D` for a single face + @unpack N, Fmask = rd + vandermonde_matrix_1D = StartUpDG.vandermonde(Line(), N, StartUpDG.nodes(Line(), N)) + rplot = LinRange(-1, 1, nvisnodes) + Vp1D = StartUpDG.vandermonde(Line(), N, rplot) / vandermonde_matrix_1D + + num_face_points = N + 1 + num_faces_total = num_faces(rd.element_type) * md.num_elements + xf, yf, uf = map(x -> reshape(view(x, Fmask, :), num_face_points, num_faces_total), + (md.xyz..., u.data)) + + num_face_plotting_points = size(Vp1D, 1) + x_mesh, y_mesh = ntuple(_ -> zeros(num_face_plotting_points, num_faces_total), 2) + u_mesh = similar(u.data, (num_face_plotting_points, num_faces_total)) + for f in 1:num_faces_total + mul!(view(x_mesh, :, f), Vp1D, view(xf, :, f)) + mul!(view(y_mesh, :, f), Vp1D, view(yf, :, f)) + mul!(view(u_mesh, :, f), Vp1D, view(uf, :, f)) + end + + return x_mesh, y_mesh, u_mesh + end + + # These methods are used internally to set the default value of the solution variables: + # - If a `cons2prim` for the given `equations` exists, use it + # - Otherwise, use `cons2cons`, which is defined for all systems of equations + digest_solution_variables(equations, solution_variables) = solution_variables + function digest_solution_variables(equations, solution_variables::Nothing) + if hasmethod(cons2prim, Tuple{AbstractVector, typeof(equations)}) + return cons2prim + else + return cons2cons + end end -end - -""" - adapt_to_mesh_level!(u_ode, semi, level) - adapt_to_mesh_level!(sol::Trixi.TrixiODESolution, level) - -Like [`adapt_to_mesh_level`](@ref), but modifies the solution and parts of the -semidiscretization (mesh and caches) in place. -""" -function adapt_to_mesh_level!(u_ode, semi, level) - # Create AMR callback with controller that refines everything towards a single level - amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first), - base_level = level) - amr_callback = AMRCallback(semi, amr_controller, interval = 0) - - # Adapt mesh until it does not change anymore - has_changed = amr_callback.affect!(u_ode, semi, 0.0, 0) - while has_changed + + """ + adapt_to_mesh_level!(u_ode, semi, level) + adapt_to_mesh_level!(sol::Trixi.TrixiODESolution, level) + + Like [`adapt_to_mesh_level`](@ref), but modifies the solution and parts of the + semidiscretization (mesh and caches) in place. + """ + function adapt_to_mesh_level!(u_ode, semi, level) + # Create AMR callback with controller that refines everything towards a single level + amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first), + base_level = level) + amr_callback = AMRCallback(semi, amr_controller, interval = 0) + + # Adapt mesh until it does not change anymore has_changed = amr_callback.affect!(u_ode, semi, 0.0, 0) - end - - return u_ode, semi -end - -function adapt_to_mesh_level!(sol::TrixiODESolution, level) - adapt_to_mesh_level!(sol.u[end], sol.prob.p, level) -end - -""" - adapt_to_mesh_level(u_ode, semi, level) - adapt_to_mesh_level(sol::Trixi.TrixiODESolution, level) - -Use the regular adaptive mesh refinement routines to adaptively refine/coarsen the solution `u_ode` -with semidiscretization `semi` towards a uniformly refined grid with refinement level `level`. The -solution and semidiscretization are copied such that the original objects remain *unaltered*. - -A convenience method accepts an ODE solution object, from which solution and semidiscretization are -extracted as needed. - -See also: [`adapt_to_mesh_level!`](@ref) -""" -function adapt_to_mesh_level(u_ode, semi, level) - # Create new semidiscretization with copy of the current mesh - mesh, _, _, _ = mesh_equations_solver_cache(semi) - new_semi = remake(semi, mesh = deepcopy(mesh)) - - return adapt_to_mesh_level!(deepcopy(u_ode), new_semi, level) -end - -function adapt_to_mesh_level(sol::TrixiODESolution, level) - adapt_to_mesh_level(sol.u[end], sol.prob.p, level) -end - -# Extract data from a 2D/3D DG solution and prepare it for visualization as a heatmap/contour plot. -# -# Returns a tuple with -# - x coordinates -# - y coordinates -# - nodal 2D data -# - x vertices for mesh visualization -# - y vertices for mesh visualization -# -# Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may -# thus be changed in future releases. -function get_data_2d(center_level_0, length_level_0, leaf_cells, coordinates, levels, - ndims, unstructured_data, n_nodes, - grid_lines = false, max_supported_level = 11, nvisnodes = nothing, - slice = :xy, point = (0.0, 0.0, 0.0)) - # Determine resolution for data interpolation - max_level = maximum(levels) - if max_level > max_supported_level - error("Maximum refinement level $max_level is higher than " * - "maximum supported level $max_supported_level") - end - max_available_nodes_per_finest_element = 2^(max_supported_level - max_level) - if nvisnodes === nothing - max_nvisnodes = 2 * n_nodes - elseif nvisnodes == 0 - max_nvisnodes = n_nodes - else - max_nvisnodes = nvisnodes - end - nvisnodes_at_max_level = min(max_available_nodes_per_finest_element, max_nvisnodes) - resolution = nvisnodes_at_max_level * 2^max_level - nvisnodes_per_level = [2^(max_level - level) * nvisnodes_at_max_level - for level in 0:max_level] - # nvisnodes_per_level is an array (accessed by "level + 1" to accommodate - # level-0-cell) that contains the number of visualization nodes for any - # refinement level to visualize on an equidistant grid - - if ndims == 3 - (unstructured_data, coordinates, levels, - center_level_0) = unstructured_3d_to_2d(unstructured_data, - coordinates, levels, length_level_0, - center_level_0, slice, - point) - end - - # Normalize element coordinates: move center to (0, 0) and domain size to [-1, 1]² - n_elements = length(levels) - normalized_coordinates = similar(coordinates) - for element_id in 1:n_elements - @views normalized_coordinates[:, element_id] .= ((coordinates[:, element_id] .- - center_level_0) ./ - (length_level_0 / 2)) - end - - # Interpolate unstructured DG data to structured data - (structured_data = unstructured2structured(unstructured_data, - normalized_coordinates, - levels, resolution, nvisnodes_per_level)) - - # Interpolate cell-centered values to node-centered values - node_centered_data = cell2node(structured_data) - - # Determine axis coordinates for contour plot - xs = collect(range(-1, 1, length = resolution + 1)) .* length_level_0 / 2 .+ - center_level_0[1] - ys = collect(range(-1, 1, length = resolution + 1)) .* length_level_0 / 2 .+ - center_level_0[2] - - # Determine element vertices to plot grid lines - if grid_lines - mesh_vertices_x, mesh_vertices_y = calc_vertices(coordinates, levels, - length_level_0) - else - mesh_vertices_x = Vector{Float64}(undef, 0) - mesh_vertices_y = Vector{Float64}(undef, 0) - end - - return xs, ys, node_centered_data, mesh_vertices_x, mesh_vertices_y -end - -# Extract data from a 1D DG solution and prepare it for visualization as a line plot. -# This returns a tuple with -# - x coordinates -# - nodal 1D data -# -# Note: This is a low-level function that is not considered as part of Trixi's interface and may -# thus be changed in future releases. -function get_data_1d(original_nodes, unstructured_data, nvisnodes) - # Get the dimensions of u; where n_vars is the number of variables, n_nodes the number of nodal values per element and n_elements the total number of elements. - n_nodes, n_elements, n_vars = size(unstructured_data) - - # Set the amount of nodes visualized according to nvisnodes. - if nvisnodes === nothing - max_nvisnodes = 2 * n_nodes - elseif nvisnodes == 0 - max_nvisnodes = n_nodes - else - @assert nvisnodes>=2 "nvisnodes must be zero or >= 2" - max_nvisnodes = nvisnodes - end - - interpolated_nodes = Array{eltype(original_nodes), 2}(undef, max_nvisnodes, - n_elements) - interpolated_data = Array{eltype(unstructured_data), 3}(undef, max_nvisnodes, - n_elements, n_vars) - - for j in 1:n_elements - # Interpolate on an equidistant grid. - interpolated_nodes[:, j] .= range(original_nodes[1, 1, j], - original_nodes[1, end, j], - length = max_nvisnodes) - end - - nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes) - nodes_out = collect(range(-1, 1, length = max_nvisnodes)) - - # Calculate vandermonde matrix for interpolation. - vandermonde = polynomial_interpolation_matrix(nodes_in, nodes_out) - - # Iterate over all variables. - for v in 1:n_vars - # Interpolate data for each element. - for element in 1:n_elements - multiply_scalar_dimensionwise!(@view(interpolated_data[:, element, v]), - vandermonde, - @view(unstructured_data[:, element, v])) + while has_changed + has_changed = amr_callback.affect!(u_ode, semi, 0.0, 0) end - end - # Return results after data is reshaped - return vec(interpolated_nodes), reshape(interpolated_data, :, n_vars), - vcat(original_nodes[1, 1, :], original_nodes[1, end, end]) -end - -# Change order of dimensions (variables are now last) and convert data to `solution_variables` -# -# Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may -# thus be changed in future releases. -function get_unstructured_data(u, solution_variables, mesh, equations, solver, cache) - if solution_variables === cons2cons - raw_data = u - n_vars = size(raw_data, 1) - else - # FIXME: Remove this comment once the implementation following it has been verified - # Reinterpret the solution array as an array of conservative variables, - # compute the solution variables via broadcasting, and reinterpret the - # result as a plain array of floating point numbers - # raw_data = Array(reinterpret(eltype(u), - # solution_variables.(reinterpret(SVector{nvariables(equations),eltype(u)}, u), - # Ref(equations)))) - # n_vars = size(raw_data, 1) - n_vars_in = nvariables(equations) - n_vars = length(solution_variables(get_node_vars(u, equations, solver), - equations)) - raw_data = Array{eltype(u)}(undef, n_vars, Base.tail(size(u))...) - reshaped_u = reshape(u, n_vars_in, :) - reshaped_r = reshape(raw_data, n_vars, :) - for idx in axes(reshaped_u, 2) - reshaped_r[:, idx] = solution_variables(get_node_vars(reshaped_u, equations, - solver, idx), - equations) + + return u_ode, semi + end + + function adapt_to_mesh_level!(sol::TrixiODESolution, level) + adapt_to_mesh_level!(sol.u[end], sol.prob.p, level) + end + + """ + adapt_to_mesh_level(u_ode, semi, level) + adapt_to_mesh_level(sol::Trixi.TrixiODESolution, level) + + Use the regular adaptive mesh refinement routines to adaptively refine/coarsen the solution `u_ode` + with semidiscretization `semi` towards a uniformly refined grid with refinement level `level`. The + solution and semidiscretization are copied such that the original objects remain *unaltered*. + + A convenience method accepts an ODE solution object, from which solution and semidiscretization are + extracted as needed. + + See also: [`adapt_to_mesh_level!`](@ref) + """ + function adapt_to_mesh_level(u_ode, semi, level) + # Create new semidiscretization with copy of the current mesh + mesh, _, _, _ = mesh_equations_solver_cache(semi) + new_semi = remake(semi, mesh = deepcopy(mesh)) + + return adapt_to_mesh_level!(deepcopy(u_ode), new_semi, level) + end + + function adapt_to_mesh_level(sol::TrixiODESolution, level) + adapt_to_mesh_level(sol.u[end], sol.prob.p, level) + end + + # Extract data from a 3D DG solution and prepare it for visualization as a heatmap/contour plot. + # + # Returns a tuple with + # - x coordinates + # - y coordinates + # - z coordinates + # - nodal 3D data + # - x vertices for mesh visualization + # - y vertices for mesh visualization + # - z vertices for mesh visualization + # + # Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may + # thus be changed in future releases. + function get_data_3d(center_level_0, length_level_0, coordinates, levels, + unstructured_data, n_nodes, + grid_lines = false, max_supported_level = 11, nvisnodes = nothing) + # Determine resolution for data interpolation + max_level = maximum(levels) + if max_level > max_supported_level + error("Maximum refinement level $max_level is higher than " * + "maximum supported level $max_supported_level") end - end - - unstructured_data = Array{eltype(raw_data)}(undef, - ntuple((d) -> nnodes(solver), - ndims(equations))..., - nelements(solver, cache), n_vars) - for variable in 1:n_vars - @views unstructured_data[.., :, variable] .= raw_data[variable, .., :] - end - - return unstructured_data -end - -# Convert cell-centered values to node-centered values by averaging over all -# four neighbors and making use of the periodicity of the solution -# -# Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may -# thus be changed in future releases. -function cell2node(cell_centered_data) - # Create temporary data structure to make the averaging algorithm as simple - # as possible (by using a ghost layer) - tmp = similar(first(cell_centered_data), size(first(cell_centered_data)) .+ (2, 2)) - - # Create output data structure - resolution_in, _ = size(first(cell_centered_data)) - resolution_out = resolution_in + 1 - node_centered_data = [Matrix{Float64}(undef, resolution_out, resolution_out) - for _ in eachindex(cell_centered_data)] - - for (cell_data, node_data) in zip(cell_centered_data, node_centered_data) - # Fill center with original data - tmp[2:(end - 1), 2:(end - 1)] .= cell_data - - # Fill sides with opposite data (periodic domain) - # x-direction - tmp[1, 2:(end - 1)] .= cell_data[end, :] - tmp[end, 2:(end - 1)] .= cell_data[1, :] - # y-direction - tmp[2:(end - 1), 1] .= cell_data[:, end] - tmp[2:(end - 1), end] .= cell_data[:, 1] - # Corners - tmp[1, 1] = cell_data[end, end] - tmp[end, 1] = cell_data[1, end] - tmp[1, end] = cell_data[end, 1] - tmp[end, end] = cell_data[1, 1] - - # Obtain node-centered value by averaging over neighboring cell-centered values - for j in 1:resolution_out - for i in 1:resolution_out - node_data[i, j] = (tmp[i, j] + - tmp[i + 1, j] + - tmp[i, j + 1] + - tmp[i + 1, j + 1]) / 4 + max_available_nodes_per_finest_element = 2^(max_supported_level - max_level) + if nvisnodes === nothing + max_nvisnodes = 2 * n_nodes + elseif nvisnodes == 0 + max_nvisnodes = n_nodes + else + max_nvisnodes = nvisnodes + end + nvisnodes_at_max_level = min(max_available_nodes_per_finest_element, max_nvisnodes) + resolution = nvisnodes_at_max_level * 2^max_level + nvisnodes_per_level = [2^(max_level - level) * nvisnodes_at_max_level + for level in 0:max_level] + # nvisnodes_per_level is an array (accessed by "level + 1" to accommodate + # level-0-cell) that contains the number of visualization nodes for any + # refinement level to visualize on an equidistant grid + + # Normalize element coordinates: move center to (0, 0) and domain size to [-1, 1]² + n_elements = length(levels) + normalized_coordinates = similar(coordinates) + for element_id in 1:n_elements + @views normalized_coordinates[:, element_id] .= ((coordinates[:, element_id] .- + center_level_0) ./ + (length_level_0 / 2)) + end + + # Interpolate unstructured DG data to structured data + (structured_data = unstructured2structured3D(unstructured_data, + normalized_coordinates, + levels, resolution, + nvisnodes_per_level)) + + # Interpolate cell-centered values to node-centered values + node_centered_data = cell2node3D(structured_data) + + # Determine axis coordinates for contour plot + xs = collect(range(-1, 1, length = resolution + 1)) .* length_level_0 / 2 .+ + center_level_0[1] + ys = collect(range(-1, 1, length = resolution + 1)) .* length_level_0 / 2 .+ + center_level_0[2] + zs = collect(range(-1, 1, length = resolution + 1)) .* length_level_0 / 2 .+ + center_level_0[3] + + # Determine element vertices to plot grid lines + if grid_lines + mesh_vertices_x, mesh_vertices_y, mesh_vertices_z = calc_vertices3D(coordinates, + levels, + length_level_0) + else + mesh_vertices_x = Vector{Float64}(undef, 0) + mesh_vertices_y = Vector{Float64}(undef, 0) + mesh_vertices_z = Vector{Float64}(undef, 0) + end + + return xs, ys, zs, node_centered_data, mesh_vertices_x, mesh_vertices_y, + mesh_vertices_z + end + + # Extract data from a 2D/3D DG solution and prepare it for visualization as a heatmap/contour plot. + # + # Returns a tuple with + # - x coordinates + # - y coordinates + # - nodal 2D data + # - x vertices for mesh visualization + # - y vertices for mesh visualization + # + # Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may + # thus be changed in future releases. + function get_data_2d(center_level_0, length_level_0, leaf_cells, coordinates, levels, + ndims, unstructured_data, n_nodes, + grid_lines = false, max_supported_level = 11, nvisnodes = nothing, + slice = :xy, point = (0.0, 0.0, 0.0)) + # Determine resolution for data interpolation + max_level = maximum(levels) + if max_level > max_supported_level + error("Maximum refinement level $max_level is higher than " * + "maximum supported level $max_supported_level") + end + max_available_nodes_per_finest_element = 2^(max_supported_level - max_level) + if nvisnodes === nothing + max_nvisnodes = 2 * n_nodes + elseif nvisnodes == 0 + max_nvisnodes = n_nodes + else + max_nvisnodes = nvisnodes + end + nvisnodes_at_max_level = min(max_available_nodes_per_finest_element, max_nvisnodes) + resolution = nvisnodes_at_max_level * 2^max_level + nvisnodes_per_level = [2^(max_level - level) * nvisnodes_at_max_level + for level in 0:max_level] + # nvisnodes_per_level is an array (accessed by "level + 1" to accommodate + # level-0-cell) that contains the number of visualization nodes for any + # refinement level to visualize on an equidistant grid + + if ndims == 3 + (unstructured_data, coordinates, levels, + center_level_0) = unstructured_3d_to_2d(unstructured_data, + coordinates, levels, length_level_0, + center_level_0, slice, + point) + end + + # Normalize element coordinates: move center to (0, 0) and domain size to [-1, 1]² + n_elements = length(levels) + normalized_coordinates = similar(coordinates) + for element_id in 1:n_elements + @views normalized_coordinates[:, element_id] .= ((coordinates[:, element_id] .- + center_level_0) ./ + (length_level_0 / 2)) + end + + # Interpolate unstructured DG data to structured data + (structured_data = unstructured2structured(unstructured_data, + normalized_coordinates, + levels, resolution, nvisnodes_per_level)) + + # Interpolate cell-centered values to node-centered values + node_centered_data = cell2node(structured_data) + + # Determine axis coordinates for contour plot + xs = collect(range(-1, 1, length = resolution + 1)) .* length_level_0 / 2 .+ + center_level_0[1] + ys = collect(range(-1, 1, length = resolution + 1)) .* length_level_0 / 2 .+ + center_level_0[2] + + # Determine element vertices to plot grid lines + if grid_lines + mesh_vertices_x, mesh_vertices_y = calc_vertices(coordinates, levels, + length_level_0) + else + mesh_vertices_x = Vector{Float64}(undef, 0) + mesh_vertices_y = Vector{Float64}(undef, 0) + end + + return xs, ys, node_centered_data, mesh_vertices_x, mesh_vertices_y + end + + # Extract data from a 1D DG solution and prepare it for visualization as a line plot. + # This returns a tuple with + # - x coordinates + # - nodal 1D data + # + # Note: This is a low-level function that is not considered as part of Trixi's interface and may + # thus be changed in future releases. + function get_data_1d(original_nodes, unstructured_data, nvisnodes) + # Get the dimensions of u; where n_vars is the number of variables, n_nodes the number of nodal values per element and n_elements the total number of elements. + n_nodes, n_elements, n_vars = size(unstructured_data) + + # Set the amount of nodes visualized according to nvisnodes. + if nvisnodes === nothing + max_nvisnodes = 2 * n_nodes + elseif nvisnodes == 0 + max_nvisnodes = n_nodes + else + @assert nvisnodes>=2 "nvisnodes must be zero or >= 2" + max_nvisnodes = nvisnodes + end + + interpolated_nodes = Array{eltype(original_nodes), 2}(undef, max_nvisnodes, + n_elements) + interpolated_data = Array{eltype(unstructured_data), 3}(undef, max_nvisnodes, + n_elements, n_vars) + + for j in 1:n_elements + # Interpolate on an equidistant grid. + interpolated_nodes[:, j] .= range(original_nodes[1, 1, j], + original_nodes[1, end, j], + length = max_nvisnodes) + end + + nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes) + nodes_out = collect(range(-1, 1, length = max_nvisnodes)) + + # Calculate vandermonde matrix for interpolation. + vandermonde = polynomial_interpolation_matrix(nodes_in, nodes_out) + + # Iterate over all variables. + for v in 1:n_vars + # Interpolate data for each element. + for element in 1:n_elements + multiply_scalar_dimensionwise!(@view(interpolated_data[:, element, v]), + vandermonde, + @view(unstructured_data[:, element, v])) end end - end - - # Transpose - for (index, data) in enumerate(node_centered_data) - node_centered_data[index] = permutedims(data) - end - - return node_centered_data -end - -# Convert 3d unstructured data to 2d data. -# Additional to the new unstructured data updated coordinates, levels and -# center coordinates are returned. -# -# Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may -# thus be changed in future releases. -function unstructured_3d_to_2d(unstructured_data, coordinates, levels, - length_level_0, center_level_0, slice, - point) - if slice === :yz - slice_dimension = 1 - other_dimensions = [2, 3] - elseif slice === :xz - slice_dimension = 2 - other_dimensions = [1, 3] - elseif slice === :xy - slice_dimension = 3 - other_dimensions = [1, 2] - else - error("illegal dimension '$slice', supported dimensions are :yz, :xz, and :xy") - end - - # Limits of domain in slice dimension - lower_limit = center_level_0[slice_dimension] - length_level_0 / 2 - upper_limit = center_level_0[slice_dimension] + length_level_0 / 2 - - @assert length(point)>=3 "Point must be three-dimensional." - if point[slice_dimension] < lower_limit || point[slice_dimension] > upper_limit - error(string("Slice plane is outside of domain.", - " point[$slice_dimension]=$(point[slice_dimension]) must be between $lower_limit and $upper_limit")) - end - - # Extract data shape information - n_nodes_in, _, _, n_elements, n_variables = size(unstructured_data) - - # Get node coordinates for DG locations on reference element - nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in) - - # New unstructured data has one dimension less. - # The redundant element ids are removed later. - @views new_unstructured_data = similar(unstructured_data[1, ..]) - - # Declare new empty arrays to fill in new coordinates and levels - new_coordinates = Array{Float64}(undef, 2, n_elements) - new_levels = Array{eltype(levels)}(undef, n_elements) - - # Counter for new element ids - new_id = 0 - - # Save vandermonde matrices in a Dict to prevent redundant generation - vandermonde_to_2d = Dict() - - # Permute dimensions such that the slice dimension is always the - # third dimension of the array. Below we can always interpolate in the - # third dimension. - if slice === :yz - unstructured_data = permutedims(unstructured_data, [2, 3, 1, 4, 5]) - elseif slice === :xz - unstructured_data = permutedims(unstructured_data, [1, 3, 2, 4, 5]) - end - - for element_id in 1:n_elements - # Distance from center to border of this element (half the length) - element_length = length_level_0 / 2^levels[element_id] - min_coordinate = coordinates[:, element_id] .- element_length / 2 - max_coordinate = coordinates[:, element_id] .+ element_length / 2 - - # Check if slice plane and current element intersect. - # The first check uses a "greater but not equal" to only match one cell if the - # slice plane lies between two cells. - # The second check is needed if the slice plane is at the upper border of - # the domain due to this. - if !((min_coordinate[slice_dimension] <= point[slice_dimension] && - max_coordinate[slice_dimension] > point[slice_dimension]) || - (point[slice_dimension] == upper_limit && - max_coordinate[slice_dimension] == upper_limit)) - # Continue for loop if they don't intersect - continue - end - - # This element is of interest - new_id += 1 - - # Add element to new coordinates and levels - new_coordinates[:, new_id] = coordinates[other_dimensions, element_id] - new_levels[new_id] = levels[element_id] - - # Construct vandermonde matrix (or load from Dict if possible) - normalized_intercept = (point[slice_dimension] - - min_coordinate[slice_dimension]) / - element_length * 2 - 1 - - if haskey(vandermonde_to_2d, normalized_intercept) - vandermonde = vandermonde_to_2d[normalized_intercept] + # Return results after data is reshaped + return vec(interpolated_nodes), reshape(interpolated_data, :, n_vars), + vcat(original_nodes[1, 1, :], original_nodes[1, end, end]) + end + + # Change order of dimensions (variables are now last) and convert data to `solution_variables` + # + # Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may + # thus be changed in future releases. + function get_unstructured_data(u, solution_variables, mesh, equations, solver, cache) + if solution_variables === cons2cons + raw_data = u + n_vars = size(raw_data, 1) else - # Generate vandermonde matrix to interpolate values at nodes_in to one value - vandermonde = polynomial_interpolation_matrix(nodes_in, - [normalized_intercept]) - vandermonde_to_2d[normalized_intercept] = vandermonde - end - - # 1D interpolation to specified slice plane - # We permuted the dimensions above such that now the dimension in which - # we will interpolate is always the third one. - for i in 1:n_nodes_in - for ii in 1:n_nodes_in - # Interpolate in the third dimension - data = unstructured_data[i, ii, :, element_id, :] - - value = multiply_dimensionwise(vandermonde, permutedims(data)) - new_unstructured_data[i, ii, new_id, :] = value[:, 1] + # FIXME: Remove this comment once the implementation following it has been verified + # Reinterpret the solution array as an array of conservative variables, + # compute the solution variables via broadcasting, and reinterpret the + # result as a plain array of floating point numbers + # raw_data = Array(reinterpret(eltype(u), + # solution_variables.(reinterpret(SVector{nvariables(equations),eltype(u)}, u), + # Ref(equations)))) + # n_vars = size(raw_data, 1) + n_vars_in = nvariables(equations) + n_vars = length(solution_variables(get_node_vars(u, equations, solver), + equations)) + raw_data = Array{eltype(u)}(undef, n_vars, Base.tail(size(u))...) + reshaped_u = reshape(u, n_vars_in, :) + reshaped_r = reshape(raw_data, n_vars, :) + for idx in axes(reshaped_u, 2) + reshaped_r[:, idx] = solution_variables(get_node_vars(reshaped_u, equations, + solver, idx), + equations) end end - end - - # Remove redundant element ids - unstructured_data = new_unstructured_data[:, :, 1:new_id, :] - new_coordinates = new_coordinates[:, 1:new_id] - new_levels = new_levels[1:new_id] - - center_level_0 = center_level_0[other_dimensions] - - return unstructured_data, new_coordinates, new_levels, center_level_0 -end - -# Convert 2d unstructured data to 1d slice and interpolate them. -function unstructured_2d_to_1d(original_nodes, unstructured_data, nvisnodes, slice, - point) - if slice === :x - slice_dimension = 2 - other_dimension = 1 - elseif slice === :y - slice_dimension = 1 - other_dimension = 2 - else - error("illegal dimension '$slice', supported dimensions are :x and :y") - end - - # Set up data structures to store new 1D data. - @views new_unstructured_data = similar(unstructured_data[1, ..]) - @views new_nodes = similar(original_nodes[1, 1, ..]) - - n_nodes_in, _, n_elements, n_variables = size(unstructured_data) - nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in) - - # Test if point lies in the domain. - lower_limit = original_nodes[1, 1, 1, 1] - upper_limit = original_nodes[1, n_nodes_in, n_nodes_in, n_elements] - - @assert length(point)>=2 "Point must be two-dimensional." - if point[slice_dimension] < lower_limit || point[slice_dimension] > upper_limit - error(string("Slice axis is outside of domain. ", - " point[$slice_dimension]=$(point[slice_dimension]) must be between $lower_limit and $upper_limit")) - end - - # Count the amount of new elements. - new_id = 0 - - # Permute dimensions so that the slice dimension is always in the correct place for later use. - if slice === :y - original_nodes = permutedims(original_nodes, [1, 3, 2, 4]) - unstructured_data = permutedims(unstructured_data, [2, 1, 3, 4]) - end - - # Iterate over all elements to find the ones that lie on the slice axis. - for element_id in 1:n_elements - min_coordinate = original_nodes[:, 1, 1, element_id] - max_coordinate = original_nodes[:, n_nodes_in, n_nodes_in, element_id] - element_length = max_coordinate - min_coordinate - - # Test if the element is on the slice axis. If not just continue with the next element. - if !((min_coordinate[slice_dimension] <= point[slice_dimension] && - max_coordinate[slice_dimension] > point[slice_dimension]) || - (point[slice_dimension] == upper_limit && - max_coordinate[slice_dimension] == upper_limit)) - continue - end - - new_id += 1 - - # Construct vandermonde matrix for interpolation of each 2D element to a 1D element. - normalized_intercept = (point[slice_dimension] - - min_coordinate[slice_dimension]) / - element_length[1] * 2 - 1 - vandermonde = polynomial_interpolation_matrix(nodes_in, normalized_intercept) - - # Interpolate to each node of new 1D element. - for v in 1:n_variables - for node in 1:n_nodes_in - new_unstructured_data[node, new_id, v] = (vandermonde * unstructured_data[node, + + unstructured_data = Array{eltype(raw_data)}(undef, + ntuple((d) -> nnodes(solver), + ndims(equations))..., + nelements(solver, cache), n_vars) + for variable in 1:n_vars + @views unstructured_data[.., :, variable] .= raw_data[variable, .., :] + end + + return unstructured_data + end + + # Convert cell-centered values to node-centered values by averaging over all + # four neighbors and making use of the periodicity of the solution + # + # Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may + # thus be changed in future releases. + function cell2node3D(cell_centered_data) + # Create temporary data structure to make the averaging algorithm as simple + # as possible (by using a ghost layer) + tmp = similar(first(cell_centered_data), + size(first(cell_centered_data)) .+ (2, 2, 2)) + + # Create output data structure + resolution_in, _, _ = size(first(cell_centered_data)) + resolution_out = resolution_in + 1 + node_centered_data = [Array{Float64, 3}(undef, resolution_out, resolution_out, + resolution_out) + for _ in eachindex(cell_centered_data)] + + for (cell_data, node_data) in zip(cell_centered_data, node_centered_data) + # Fill center with original data + tmp[2:(end - 1), 2:(end - 1), 2:(end - 1)] .= cell_data + + # Fill faces with opposite data (periodic domain) + # x-direction + tmp[1, 2:(end - 1), 2:(end - 1)] .= cell_data[end, :, :] + tmp[end, 2:(end - 1), 2:(end - 1)] .= cell_data[1, :, :] + # y-direction + tmp[2:(end - 1), 1, 2:(end - 1)] .= cell_data[:, end, :] + tmp[2:(end - 1), end, 2:(end - 1)] .= cell_data[:, 1, :] + # z-direction + tmp[2:(end - 1), 2:(end - 1), 1] .= cell_data[:, :, end] + tmp[2:(end - 1), 2:(end - 1), end] .= cell_data[:, :, 1] + + # Edges + # x-direction + tmp[2:(end - 1), 1, 1] .= cell_data[:, end, 1] + tmp[2:(end - 1), end, 1] .= cell_data[:, 1, 1] + tmp[2:(end - 1), 1, end] .= cell_data[:, end, end] + tmp[2:(end - 1), end, end] .= cell_data[:, 1, end] + # y-direction + tmp[1, 2:(end - 1), 1] .= cell_data[end, :, 1] + tmp[end, 2:(end - 1), 1] .= cell_data[1, :, 1] + tmp[1, 2:(end - 1), end] .= cell_data[end, :, end] + tmp[end, 2:(end - 1), end] .= cell_data[1, :, end] + # z-direction + tmp[1, 1, 2:(end - 1)] .= cell_data[end, 1, :] + tmp[end, 1, 2:(end - 1)] .= cell_data[1, 1, :] + tmp[1, end, 2:(end - 1)] .= cell_data[end, end, :] + tmp[end, end, 2:(end - 1)] .= cell_data[1, end, :] + + # Corners + tmp[1, 1, 1] = cell_data[end, end, end] + tmp[end, 1, 1] = cell_data[1, end, end] + tmp[1, end, 1] = cell_data[end, 1, end] + tmp[end, end, 1] = cell_data[1, 1, end] + tmp[1, 1, end] = cell_data[end, end, 1] + tmp[end, 1, end] = cell_data[1, end, 1] + tmp[1, end, end] = cell_data[end, 1, 1] + tmp[end, end, end] = cell_data[1, 1, 1] + + # Obtain node-centered value by averaging over neighboring cell-centered values + for k in 1:resolution_out + for j in 1:resolution_out + for i in 1:resolution_out + node_data[i, j, k] = (tmp[i, j, k] + + tmp[i + 1, j, k] + + tmp[i, j + 1, k] + + tmp[i + 1, j + 1, k] + + tmp[i, j, k + 1] + + tmp[i + 1, j, k + 1] + + tmp[i, j + 1, k + 1] + + tmp[i + 1, j + 1, k + 1]) / 8 + end + end + end + end + + # Transpose + # TODO ??? + for (index, data) in enumerate(node_centered_data) + node_centered_data[index] = permutedims(data, (3, 2, 1)) + end + + return node_centered_data + end + + # Convert cell-centered values to node-centered values by averaging over all + # four neighbors and making use of the periodicity of the solution + # + # Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may + # thus be changed in future releases. + function cell2node(cell_centered_data) + # Create temporary data structure to make the averaging algorithm as simple + # as possible (by using a ghost layer) + tmp = similar(first(cell_centered_data), size(first(cell_centered_data)) .+ (2, 2)) + + # Create output data structure + resolution_in, _ = size(first(cell_centered_data)) + resolution_out = resolution_in + 1 + node_centered_data = [Matrix{Float64}(undef, resolution_out, resolution_out) + for _ in eachindex(cell_centered_data)] + + for (cell_data, node_data) in zip(cell_centered_data, node_centered_data) + # Fill center with original data + tmp[2:(end - 1), 2:(end - 1)] .= cell_data + + # Fill sides with opposite data (periodic domain) + # x-direction + tmp[1, 2:(end - 1)] .= cell_data[end, :] + tmp[end, 2:(end - 1)] .= cell_data[1, :] + # y-direction + tmp[2:(end - 1), 1] .= cell_data[:, end] + tmp[2:(end - 1), end] .= cell_data[:, 1] + # Corners + tmp[1, 1] = cell_data[end, end] + tmp[end, 1] = cell_data[1, end] + tmp[1, end] = cell_data[end, 1] + tmp[end, end] = cell_data[1, 1] + + # Obtain node-centered value by averaging over neighboring cell-centered values + for j in 1:resolution_out + for i in 1:resolution_out + node_data[i, j] = (tmp[i, j] + + tmp[i + 1, j] + + tmp[i, j + 1] + + tmp[i + 1, j + 1]) / 4 + end + end + end + + # Transpose + for (index, data) in enumerate(node_centered_data) + node_centered_data[index] = permutedims(data) + end + + return node_centered_data + end + + # Convert 3d unstructured data to 2d data. + # Additional to the new unstructured data updated coordinates, levels and + # center coordinates are returned. + # + # Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may + # thus be changed in future releases. + function unstructured_3d_to_2d(unstructured_data, coordinates, levels, + length_level_0, center_level_0, slice, + point) + if slice === :yz + slice_dimension = 1 + other_dimensions = [2, 3] + elseif slice === :xz + slice_dimension = 2 + other_dimensions = [1, 3] + elseif slice === :xy + slice_dimension = 3 + other_dimensions = [1, 2] + else + error("illegal dimension '$slice', supported dimensions are :yz, :xz, and :xy") + end + + # Limits of domain in slice dimension + lower_limit = center_level_0[slice_dimension] - length_level_0 / 2 + upper_limit = center_level_0[slice_dimension] + length_level_0 / 2 + + @assert length(point)>=3 "Point must be three-dimensional." + if point[slice_dimension] < lower_limit || point[slice_dimension] > upper_limit + error(string("Slice plane is outside of domain.", + " point[$slice_dimension]=$(point[slice_dimension]) must be between $lower_limit and $upper_limit")) + end + + # Extract data shape information + n_nodes_in, _, _, n_elements, n_variables = size(unstructured_data) + + # Get node coordinates for DG locations on reference element + nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in) + + # New unstructured data has one dimension less. + # The redundant element ids are removed later. + @views new_unstructured_data = similar(unstructured_data[1, ..]) + + # Declare new empty arrays to fill in new coordinates and levels + new_coordinates = Array{Float64}(undef, 2, n_elements) + new_levels = Array{eltype(levels)}(undef, n_elements) + + # Counter for new element ids + new_id = 0 + + # Save vandermonde matrices in a Dict to prevent redundant generation + vandermonde_to_2d = Dict() + + # Permute dimensions such that the slice dimension is always the + # third dimension of the array. Below we can always interpolate in the + # third dimension. + if slice === :yz + unstructured_data = permutedims(unstructured_data, [2, 3, 1, 4, 5]) + elseif slice === :xz + unstructured_data = permutedims(unstructured_data, [1, 3, 2, 4, 5]) + end + + for element_id in 1:n_elements + # Distance from center to border of this element (half the length) + element_length = length_level_0 / 2^levels[element_id] + min_coordinate = coordinates[:, element_id] .- element_length / 2 + max_coordinate = coordinates[:, element_id] .+ element_length / 2 + + # Check if slice plane and current element intersect. + # The first check uses a "greater but not equal" to only match one cell if the + # slice plane lies between two cells. + # The second check is needed if the slice plane is at the upper border of + # the domain due to this. + if !((min_coordinate[slice_dimension] <= point[slice_dimension] && + max_coordinate[slice_dimension] > point[slice_dimension]) || + (point[slice_dimension] == upper_limit && + max_coordinate[slice_dimension] == upper_limit)) + # Continue for loop if they don't intersect + continue + end + + # This element is of interest + new_id += 1 + + # Add element to new coordinates and levels + new_coordinates[:, new_id] = coordinates[other_dimensions, element_id] + new_levels[new_id] = levels[element_id] + + # Construct vandermonde matrix (or load from Dict if possible) + normalized_intercept = (point[slice_dimension] - + min_coordinate[slice_dimension]) / + element_length * 2 - 1 + + if haskey(vandermonde_to_2d, normalized_intercept) + vandermonde = vandermonde_to_2d[normalized_intercept] + else + # Generate vandermonde matrix to interpolate values at nodes_in to one value + vandermonde = polynomial_interpolation_matrix(nodes_in, + [normalized_intercept]) + vandermonde_to_2d[normalized_intercept] = vandermonde + end + + # 1D interpolation to specified slice plane + # We permuted the dimensions above such that now the dimension in which + # we will interpolate is always the third one. + for i in 1:n_nodes_in + for ii in 1:n_nodes_in + # Interpolate in the third dimension + data = unstructured_data[i, ii, :, element_id, :] + + value = multiply_dimensionwise(vandermonde, permutedims(data)) + new_unstructured_data[i, ii, new_id, :] = value[:, 1] + end + end + end + + # Remove redundant element ids + unstructured_data = new_unstructured_data[:, :, 1:new_id, :] + new_coordinates = new_coordinates[:, 1:new_id] + new_levels = new_levels[1:new_id] + + center_level_0 = center_level_0[other_dimensions] + + return unstructured_data, new_coordinates, new_levels, center_level_0 + end + + # Convert 2d unstructured data to 1d slice and interpolate them. + function unstructured_2d_to_1d(original_nodes, unstructured_data, nvisnodes, slice, + point) + if slice === :x + slice_dimension = 2 + other_dimension = 1 + elseif slice === :y + slice_dimension = 1 + other_dimension = 2 + else + error("illegal dimension '$slice', supported dimensions are :x and :y") + end + + # Set up data structures to store new 1D data. + @views new_unstructured_data = similar(unstructured_data[1, ..]) + @views new_nodes = similar(original_nodes[1, 1, ..]) + + n_nodes_in, _, n_elements, n_variables = size(unstructured_data) + nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in) + + # Test if point lies in the domain. + lower_limit = original_nodes[1, 1, 1, 1] + upper_limit = original_nodes[1, n_nodes_in, n_nodes_in, n_elements] + + @assert length(point)>=2 "Point must be two-dimensional." + if point[slice_dimension] < lower_limit || point[slice_dimension] > upper_limit + error(string("Slice axis is outside of domain. ", + " point[$slice_dimension]=$(point[slice_dimension]) must be between $lower_limit and $upper_limit")) + end + + # Count the amount of new elements. + new_id = 0 + + # Permute dimensions so that the slice dimension is always in the correct place for later use. + if slice === :y + original_nodes = permutedims(original_nodes, [1, 3, 2, 4]) + unstructured_data = permutedims(unstructured_data, [2, 1, 3, 4]) + end + + # Iterate over all elements to find the ones that lie on the slice axis. + for element_id in 1:n_elements + min_coordinate = original_nodes[:, 1, 1, element_id] + max_coordinate = original_nodes[:, n_nodes_in, n_nodes_in, element_id] + element_length = max_coordinate - min_coordinate + + # Test if the element is on the slice axis. If not just continue with the next element. + if !((min_coordinate[slice_dimension] <= point[slice_dimension] && + max_coordinate[slice_dimension] > point[slice_dimension]) || + (point[slice_dimension] == upper_limit && + max_coordinate[slice_dimension] == upper_limit)) + continue + end + + new_id += 1 + + # Construct vandermonde matrix for interpolation of each 2D element to a 1D element. + normalized_intercept = (point[slice_dimension] - + min_coordinate[slice_dimension]) / + element_length[1] * 2 - 1 + vandermonde = polynomial_interpolation_matrix(nodes_in, normalized_intercept) + + # Interpolate to each node of new 1D element. + for v in 1:n_variables + for node in 1:n_nodes_in + new_unstructured_data[node, new_id, v] = (vandermonde * unstructured_data[node, + :, + element_id, + v])[1] + end + end + + new_nodes[:, new_id] = original_nodes[other_dimension, :, 1, element_id] + end + + return get_data_1d(reshape(new_nodes[:, 1:new_id], 1, n_nodes_in, new_id), + new_unstructured_data[:, 1:new_id, :], nvisnodes) + end + + # Calculate the arc length of a curve given by ndims x npoints point coordinates (piece-wise linear approximation) + function calc_arc_length(coordinates) + n_points = size(coordinates)[2] + arc_length = zeros(n_points) + for i in 1:(n_points - 1) + arc_length[i + 1] = arc_length[i] + + sqrt(sum((coordinates[:, i] - coordinates[:, i + 1]) .^ 2)) + end + return arc_length + end + + # Convert 2d unstructured data to 1d data at given curve. + function unstructured_2d_to_1d_curve(original_nodes, unstructured_data, nvisnodes, + curve, mesh, solver, cache) + n_points_curve = size(curve)[2] + n_nodes, _, n_elements, n_variables = size(unstructured_data) + nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes) + + # Check if input is correct. + min = original_nodes[:, 1, 1, 1] + max = max_coordinate = original_nodes[:, n_nodes, n_nodes, n_elements] + @assert size(curve)==(2, size(curve)[2]) "Coordinates along curve must be 2xn dimensional." + for element in 1:n_points_curve + @assert (prod(vcat(curve[:, n_points_curve] .>= min, + curve[:, n_points_curve] + .<= + max))) "Some coordinates from `curve` are outside of the domain.." + end + + # Set nodes according to the length of the curve. + arc_length = calc_arc_length(curve) + + # Setup data structures. + data_on_curve = Array{Float64}(undef, n_points_curve, n_variables) + temp_data = Array{Float64}(undef, n_nodes, n_points_curve, n_variables) + + # For each coordinate find the corresponding element with its id. + element_ids = get_elements_by_coordinates(curve, mesh, solver, cache) + + # Iterate over all found elements. + for element in 1:n_points_curve + min_coordinate = original_nodes[:, 1, 1, element_ids[element]] + max_coordinate = original_nodes[:, n_nodes, n_nodes, element_ids[element]] + element_length = max_coordinate - min_coordinate + + normalized_coordinates = (curve[:, element] - min_coordinate) / + element_length[1] * 2 .- 1 + + # Interpolate to a single point in each element. + vandermonde_x = polynomial_interpolation_matrix(nodes_in, + normalized_coordinates[1]) + vandermonde_y = polynomial_interpolation_matrix(nodes_in, + normalized_coordinates[2]) + for v in 1:n_variables + for i in 1:n_nodes + temp_data[i, element, v] = (vandermonde_y * unstructured_data[i, :, + element_ids[element], + v])[1] + end + data_on_curve[element, v] = (vandermonde_x * temp_data[:, element, v])[] + end + end + + return arc_length, data_on_curve, nothing + end + + # Convert a PlotData2DTriangulate object to a 1d data along given curve. + function unstructured_2d_to_1d_curve(pd, input_curve, slice, point, nvisnodes) + + # If no curve is defined, create a axis curve. + if input_curve === nothing + input_curve = axis_curve(pd.x, pd.y, nothing, slice, point, nvisnodes) + end + + @assert size(input_curve, 1)==2 "Input 'curve' must be 2xn dimensional." + + # For each coordinate find the corresponding triangle with its ids. + ids_by_coordinates = get_ids_by_coordinates(input_curve, pd) + found_coordinates = ids_by_coordinates[:, 1] .!= nothing + + @assert found_coordinates!=zeros(size(input_curve, 2)) "No points of 'curve' are inside of the solutions domain." + + # These hold the ids of the elements and triangles the points of the curve sit in. + element_ids = @view ids_by_coordinates[found_coordinates, 1] + triangle_ids = @view ids_by_coordinates[found_coordinates, 2] + + # Shorten the curve, so that it contains only point that were found. + curve = @view input_curve[:, found_coordinates] + + n_variables = length(pd.data[1, 1]) + n_points_curve = size(curve, 2) + + # Set nodes according to the length of the curve. + arc_length = calc_arc_length(curve) + + # Setup data structures. + data_on_curve = Array{Float64}(undef, n_points_curve, n_variables) + + # Iterate over all points on the curve. + for point in 1:n_points_curve + element = @view element_ids[point] + triangle = @view pd.t[triangle_ids[point], :] + for v in 1:n_variables + # Get the x and y coordinates of the corners of given triangle. + x_coordinates_triangle = SVector{3}(pd.x[triangle, element]) + y_coordinates_triangle = SVector{3}(pd.y[triangle, element]) + + # Extract solutions values in corners of the triangle. + values_triangle = SVector{3}(getindex.(view(pd.data, triangle, element), v)) + + # Linear interpolation in each triangle to the points on the curve. + data_on_curve[point, v] = triangle_interpolation(x_coordinates_triangle, + y_coordinates_triangle, + values_triangle, + curve[:, point]) + end + end + + return arc_length, data_on_curve, nothing + end + + # Convert 3d unstructured data to 1d data at given curve. + function unstructured_3d_to_1d_curve(original_nodes, unstructured_data, nvisnodes, + curve, mesh, solver, cache) + n_points_curve = size(curve)[2] + n_nodes, _, _, n_elements, n_variables = size(unstructured_data) + nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes) + + # Check if input is correct. + min = original_nodes[:, 1, 1, 1, 1] + max = max_coordinate = original_nodes[:, n_nodes, n_nodes, n_nodes, n_elements] + @assert size(curve)==(3, n_points_curve) "Coordinates along curve must be 3xn dimensional." + for element in 1:n_points_curve + @assert (prod(vcat(curve[:, n_points_curve] .>= min, + curve[:, n_points_curve] + .<= + max))) "Some coordinates from `curve` are outside of the domain.." + end + + # Set nodes according to the length of the curve. + arc_length = calc_arc_length(curve) + + # Setup data structures. + data_on_curve = Array{Float64}(undef, n_points_curve, n_variables) + temp_data = Array{Float64}(undef, n_nodes, n_nodes + 1, n_points_curve, n_variables) + + # For each coordinate find the corresponding element with its id. + element_ids = get_elements_by_coordinates(curve, mesh, solver, cache) + + # Iterate over all found elements. + for element in 1:n_points_curve + min_coordinate = original_nodes[:, 1, 1, 1, element_ids[element]] + max_coordinate = original_nodes[:, n_nodes, n_nodes, n_nodes, + element_ids[element]] + element_length = max_coordinate - min_coordinate + + normalized_coordinates = (curve[:, element] - min_coordinate) / + element_length[1] * 2 .- 1 + + # Interpolate to a single point in each element. + vandermonde_x = polynomial_interpolation_matrix(nodes_in, + normalized_coordinates[1]) + vandermonde_y = polynomial_interpolation_matrix(nodes_in, + normalized_coordinates[2]) + vandermonde_z = polynomial_interpolation_matrix(nodes_in, + normalized_coordinates[3]) + for v in 1:n_variables + for i in 1:n_nodes + for ii in 1:n_nodes + temp_data[i, ii, element, v] = (vandermonde_z * unstructured_data[i, + ii, :, - element_id, + element_ids[element], v])[1] + end + temp_data[i, n_nodes + 1, element, v] = (vandermonde_y * temp_data[i, + 1:n_nodes, + element, + v])[1] + end + data_on_curve[element, v] = (vandermonde_x * temp_data[:, n_nodes + 1, + element, v])[1] end end - - new_nodes[:, new_id] = original_nodes[other_dimension, :, 1, element_id] - end - - return get_data_1d(reshape(new_nodes[:, 1:new_id], 1, n_nodes_in, new_id), - new_unstructured_data[:, 1:new_id, :], nvisnodes) -end - -# Calculate the arc length of a curve given by ndims x npoints point coordinates (piece-wise linear approximation) -function calc_arc_length(coordinates) - n_points = size(coordinates)[2] - arc_length = zeros(n_points) - for i in 1:(n_points - 1) - arc_length[i + 1] = arc_length[i] + - sqrt(sum((coordinates[:, i] - coordinates[:, i + 1]) .^ 2)) - end - return arc_length -end - -# Convert 2d unstructured data to 1d data at given curve. -function unstructured_2d_to_1d_curve(original_nodes, unstructured_data, nvisnodes, - curve, mesh, solver, cache) - n_points_curve = size(curve)[2] - n_nodes, _, n_elements, n_variables = size(unstructured_data) - nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes) - - # Check if input is correct. - min = original_nodes[:, 1, 1, 1] - max = max_coordinate = original_nodes[:, n_nodes, n_nodes, n_elements] - @assert size(curve)==(2, size(curve)[2]) "Coordinates along curve must be 2xn dimensional." - for element in 1:n_points_curve - @assert (prod(vcat(curve[:, n_points_curve] .>= min, - curve[:, n_points_curve] - .<= - max))) "Some coordinates from `curve` are outside of the domain.." + + return arc_length, data_on_curve, nothing + end + + # Convert 3d unstructured data from a general mesh to 1d data at given curve. + function unstructured_3d_to_1d_curve(nodes, data, curve, slice, point, nvisnodes) + # If no curve is defined, create a axis curve. + if curve === nothing + curve = axis_curve(nodes[1, :, :, :, :], nodes[2, :, :, :, :], + nodes[3, :, :, :, :], slice, point, nvisnodes) + end + + # Set up data structure. + n_points_curve = size(curve, 2) + n_variables = size(data, 1) + data_on_curve = Array{Float64}(undef, n_points_curve, n_variables) + + # Iterate over every point on the curve and determine the solutions value at given point. + for i in 1:n_points_curve + @views data_on_curve[i, :] .= get_value_at_point(curve[:, i], nodes, data) + end + + mesh_vertices_x = nothing + + return calc_arc_length(curve), data_on_curve, mesh_vertices_x + end + + # Check if the first 'amount'-many points can still form a valid tetrahedron. + function is_valid_tetrahedron(amount, coordinates; tol = 10^-4) + a = coordinates[:, 1] + b = coordinates[:, 2] + c = coordinates[:, 3] + d = coordinates[:, 4] + if amount == 2 # If two points are the same, then no tetrahedron can be formed. + return !(isapprox(a, b; atol = tol)) + elseif amount == 3 # Check if three points are on the same line. + return !on_the_same_line(a, b, c; tol = tol) + elseif amount == 4 # Check if four points form a tetrahedron. + A = hcat(coordinates[1, :], coordinates[2, :], coordinates[3, :], + SVector(1, 1, 1, 1)) + return !isapprox(det(A), 0; atol = tol) + else # With one point a tetrahedron can always be formed. + return true + end end - - # Set nodes according to the length of the curve. - arc_length = calc_arc_length(curve) - - # Setup data structures. - data_on_curve = Array{Float64}(undef, n_points_curve, n_variables) - temp_data = Array{Float64}(undef, n_nodes, n_points_curve, n_variables) - - # For each coordinate find the corresponding element with its id. - element_ids = get_elements_by_coordinates(curve, mesh, solver, cache) - - # Iterate over all found elements. - for element in 1:n_points_curve - min_coordinate = original_nodes[:, 1, 1, element_ids[element]] - max_coordinate = original_nodes[:, n_nodes, n_nodes, element_ids[element]] - element_length = max_coordinate - min_coordinate - - normalized_coordinates = (curve[:, element] - min_coordinate) / - element_length[1] * 2 .- 1 - - # Interpolate to a single point in each element. - vandermonde_x = polynomial_interpolation_matrix(nodes_in, - normalized_coordinates[1]) - vandermonde_y = polynomial_interpolation_matrix(nodes_in, - normalized_coordinates[2]) + + # Check if three given 3D-points are on the same line. + function on_the_same_line(a, b, c; tol = 10^-4) + # Calculate the intersection of the a-b-axis at x=0. + if b[1] == 0 + intersect_a_b = b + else + intersect_a_b = a - b .* (a[1] / b[1]) + end + # Calculate the intersection of the a-c-axis at x=0. + if c[1] == 0 + intersect_a_c = c + else + intersect_a_c = a - c .* (a[1] / c[1]) + end + return isapprox(intersect_a_b, intersect_a_c; atol = tol) + end + + # Interpolate from four corners of a tetrahedron to a single point. + function tetrahedron_interpolation(x_coordinates_in, y_coordinates_in, z_coordinates_in, + values_in, coordinate_out) + A = hcat(x_coordinates_in, y_coordinates_in, z_coordinates_in, SVector(1, 1, 1, 1)) + c = A \ values_in + return c[1] * coordinate_out[1] + c[2] * coordinate_out[2] + + c[3] * coordinate_out[3] + c[4] + end + + # Calculate the distances from every entry in node to given point. + function distances_from_single_point(nodes, point) + _, n_nodes, _, _, n_elements = size(nodes) + shifted_data = nodes .- point + distances = zeros(n_nodes, n_nodes, n_nodes, n_elements) + + # Iterate over every entry. + for element in 1:n_elements + for x in 1:n_nodes + for y in 1:n_nodes + for z in 1:n_nodes + distances[x, y, z, element] = norm(shifted_data[:, x, y, z, + element]) + end + end + end + end + return distances + end + + # Interpolate the data on given nodes to a single value at given point. + function get_value_at_point(point, nodes, data) + # Set up ata structures. + n_variables, n_x_nodes, n_y_nodes, n_z_nodes, _ = size(data) + distances = distances_from_single_point(nodes, point) + maximum_distance = maximum(distances) + + coordinates_tetrahedron = Array{Float64, 2}(undef, 3, 4) + value_tetrahedron = Array{Float64}(undef, n_variables, 4) + + index = argmin(distances) + + # If the point sits exactly on a node, no interpolation is needed. + if nodes[:, index[1], index[2], index[3], index[4]] == point + return data[1, index[1], index[2], index[3], index[4]] + end + + @views coordinates_tetrahedron[:, 1] = nodes[:, index[1], index[2], index[3], + index[4]] + @views value_tetrahedron[:, 1] = data[:, index[1], index[2], index[3], index[4]] + + # Restrict the interpolation to the closest element only. + closest_element = index[4] + @views element_distances = distances[:, :, :, closest_element] + + # Find a tetrahedron, which is given by four corners, to interpolate from. + for i in 1:4 + # Iterate until a valid tetrahedron is found. + while true + index = argmin(element_distances) + element_distances[index[1], index[2], index[3]] = maximum_distance + + @views coordinates_tetrahedron[:, i] = nodes[:, index[1], index[2], + index[3], closest_element] + @views value_tetrahedron[:, i] = data[:, index[1], index[2], index[3], + closest_element] + + # Look for another point if current tetrahedron is not valid. + if is_valid_tetrahedron(i, coordinates_tetrahedron) + break + end + end + end + + # Interpolate from tetrahedron to given point. + value_at_point = Array{Float64}(undef, n_variables) for v in 1:n_variables - for i in 1:n_nodes - temp_data[i, element, v] = (vandermonde_y * unstructured_data[i, :, - element_ids[element], - v])[1] + value_at_point[v] = tetrahedron_interpolation(coordinates_tetrahedron[1, :], + coordinates_tetrahedron[2, :], + coordinates_tetrahedron[3, :], + value_tetrahedron[v, :], point) + end + + return value_at_point + end + + # Convert 3d unstructured data to 1d slice and interpolate them. + function unstructured_3d_to_1d(original_nodes, unstructured_data, nvisnodes, slice, + point) + if slice === :x + slice_dimension = 1 + other_dimensions = [2, 3] + elseif slice === :y + slice_dimension = 2 + other_dimensions = [1, 3] + elseif slice === :z + slice_dimension = 3 + other_dimensions = [1, 2] + else + error("illegal dimension '$slice', supported dimensions are :x, :y and :z") + end + + # Set up data structures to store new 1D data. + @views new_unstructured_data = similar(unstructured_data[1, 1, ..]) + @views temp_unstructured_data = similar(unstructured_data[1, ..]) + @views new_nodes = similar(original_nodes[1, 1, 1, ..]) + + n_nodes_in, _, _, n_elements, n_variables = size(unstructured_data) + nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in) + + # Test if point lies in the domain. + lower_limit = original_nodes[1, 1, 1, 1, 1] + upper_limit = original_nodes[1, n_nodes_in, n_nodes_in, n_nodes_in, n_elements] + + @assert length(point)>=3 "Point must be three-dimensional." + if prod(point[other_dimensions] .< lower_limit) || + prod(point[other_dimensions] .> upper_limit) + error(string("Slice axis is outside of domain. ", + " point[$other_dimensions]=$(point[other_dimensions]) must be between $lower_limit and $upper_limit")) + end + + # Count the amount of new elements. + new_id = 0 + + # Permute dimensions so that the slice dimensions are always the in correct places for later use. + if slice === :x + original_nodes = permutedims(original_nodes, [1, 3, 4, 2, 5]) + unstructured_data = permutedims(unstructured_data, [2, 3, 1, 4, 5]) + elseif slice === :y + original_nodes = permutedims(original_nodes, [1, 2, 4, 3, 5]) + unstructured_data = permutedims(unstructured_data, [1, 3, 2, 4, 5]) + end + + # Iterate over all elements to find the ones that lie on the slice axis. + for element_id in 1:n_elements + min_coordinate = original_nodes[:, 1, 1, 1, element_id] + max_coordinate = original_nodes[:, n_nodes_in, n_nodes_in, n_nodes_in, + element_id] + element_length = max_coordinate - min_coordinate + + # Test if the element is on the slice axis. If not just continue with the next element. + if !((prod(min_coordinate[other_dimensions] .<= point[other_dimensions]) && + prod(max_coordinate[other_dimensions] .> point[other_dimensions])) || + (point[other_dimensions] == upper_limit && + prod(max_coordinate[other_dimensions] .== upper_limit))) + continue + end + + new_id += 1 + + # Construct vandermonde matrix for interpolation of each 2D element to a 1D element. + normalized_intercept = (point[other_dimensions] .- + min_coordinate[other_dimensions]) / + element_length[1] * 2 .- 1 + vandermonde_i = polynomial_interpolation_matrix(nodes_in, + normalized_intercept[1]) + vandermonde_ii = polynomial_interpolation_matrix(nodes_in, + normalized_intercept[2]) + + # Interpolate to each node of new 1D element. + for v in 1:n_variables + for i in 1:n_nodes_in + for ii in 1:n_nodes_in + temp_unstructured_data[i, ii, new_id, v] = (vandermonde_ii * unstructured_data[ii, + :, + i, + element_id, + v])[1] + end + new_unstructured_data[i, new_id, v] = (vandermonde_i * temp_unstructured_data[i, + :, + new_id, + v])[1] + end end - data_on_curve[element, v] = (vandermonde_x * temp_data[:, element, v])[] + + new_nodes[:, new_id] = original_nodes[slice_dimension, 1, 1, :, element_id] end - end - - return arc_length, data_on_curve, nothing -end - -# Convert a PlotData2DTriangulate object to a 1d data along given curve. -function unstructured_2d_to_1d_curve(pd, input_curve, slice, point, nvisnodes) - - # If no curve is defined, create a axis curve. - if input_curve === nothing - input_curve = axis_curve(pd.x, pd.y, nothing, slice, point, nvisnodes) - end - - @assert size(input_curve, 1)==2 "Input 'curve' must be 2xn dimensional." - - # For each coordinate find the corresponding triangle with its ids. - ids_by_coordinates = get_ids_by_coordinates(input_curve, pd) - found_coordinates = ids_by_coordinates[:, 1] .!= nothing - - @assert found_coordinates!=zeros(size(input_curve, 2)) "No points of 'curve' are inside of the solutions domain." - - # These hold the ids of the elements and triangles the points of the curve sit in. - element_ids = @view ids_by_coordinates[found_coordinates, 1] - triangle_ids = @view ids_by_coordinates[found_coordinates, 2] - - # Shorten the curve, so that it contains only point that were found. - curve = @view input_curve[:, found_coordinates] - - n_variables = length(pd.data[1, 1]) - n_points_curve = size(curve, 2) - - # Set nodes according to the length of the curve. - arc_length = calc_arc_length(curve) - - # Setup data structures. - data_on_curve = Array{Float64}(undef, n_points_curve, n_variables) - - # Iterate over all points on the curve. - for point in 1:n_points_curve - element = @view element_ids[point] - triangle = @view pd.t[triangle_ids[point], :] + + return get_data_1d(reshape(new_nodes[:, 1:new_id], 1, n_nodes_in, new_id), + new_unstructured_data[:, 1:new_id, :], nvisnodes) + end + + # Interpolate unstructured DG data to structured data (cell-centered) + # + # This function takes DG data in an unstructured, Cartesian layout and converts it to a uniformly + # distributed 3D layout. + # + # Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may + # thus be changed in future releases. + function unstructured2structured3D(unstructured_data, normalized_coordinates, + levels, resolution, nvisnodes_per_level) + # Extract data shape information + n_nodes_in, _, _, n_elements, n_variables = size(unstructured_data) + + # Get node coordinates for DG locations on reference element + nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in) + + # Calculate interpolation vandermonde matrices for each level + max_level = length(nvisnodes_per_level) - 1 + vandermonde_per_level = [] + for l in 0:max_level + n_nodes_out = nvisnodes_per_level[l + 1] + dx = 2 / n_nodes_out + nodes_out = collect(range(-1 + dx / 2, 1 - dx / 2, length = n_nodes_out)) + push!(vandermonde_per_level, + polynomial_interpolation_matrix(nodes_in, nodes_out)) + end + + # For each element, calculate index position at which to insert data in global data structure + lower_left_index = element2index3D(normalized_coordinates, levels, resolution, + nvisnodes_per_level) + + # Create output data structure + structured = [Array{Float64, 3}(undef, resolution, resolution, resolution) + for _ in 1:n_variables] + + # For each variable, interpolate element data and store to global data structure for v in 1:n_variables - # Get the x and y coordinates of the corners of given triangle. - x_coordinates_triangle = SVector{3}(pd.x[triangle, element]) - y_coordinates_triangle = SVector{3}(pd.y[triangle, element]) - - # Extract solutions values in corners of the triangle. - values_triangle = SVector{3}(getindex.(view(pd.data, triangle, element), v)) - - # Linear interpolation in each triangle to the points on the curve. - data_on_curve[point, v] = triangle_interpolation(x_coordinates_triangle, - y_coordinates_triangle, - values_triangle, - curve[:, point]) + # Reshape data array for use in multiply_dimensionwise function + reshaped_data = reshape(unstructured_data[:, :, :, :, v], 1, n_nodes_in, + n_nodes_in, n_nodes_in, n_elements) + + for element_id in 1:n_elements + # Extract level for convenience + level = levels[element_id] + + # Determine target indices + n_nodes_out = nvisnodes_per_level[level + 1] + first = lower_left_index[:, element_id] + last = first .+ (n_nodes_out - 1) + + # Interpolate data + vandermonde = vandermonde_per_level[level + 1] + structured[v][first[1]:last[1], first[2]:last[2], first[3]:last[3]] .= (reshape(multiply_dimensionwise(vandermonde, + reshaped_data[:, + :, + :, + :, + element_id]), + n_nodes_out, + n_nodes_out, + n_nodes_out)) + end end - end - - return arc_length, data_on_curve, nothing -end - -# Convert 3d unstructured data to 1d data at given curve. -function unstructured_3d_to_1d_curve(original_nodes, unstructured_data, nvisnodes, - curve, mesh, solver, cache) - n_points_curve = size(curve)[2] - n_nodes, _, _, n_elements, n_variables = size(unstructured_data) - nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes) - - # Check if input is correct. - min = original_nodes[:, 1, 1, 1, 1] - max = max_coordinate = original_nodes[:, n_nodes, n_nodes, n_nodes, n_elements] - @assert size(curve)==(3, n_points_curve) "Coordinates along curve must be 3xn dimensional." - for element in 1:n_points_curve - @assert (prod(vcat(curve[:, n_points_curve] .>= min, - curve[:, n_points_curve] - .<= - max))) "Some coordinates from `curve` are outside of the domain.." - end - - # Set nodes according to the length of the curve. - arc_length = calc_arc_length(curve) - - # Setup data structures. - data_on_curve = Array{Float64}(undef, n_points_curve, n_variables) - temp_data = Array{Float64}(undef, n_nodes, n_nodes + 1, n_points_curve, n_variables) - - # For each coordinate find the corresponding element with its id. - element_ids = get_elements_by_coordinates(curve, mesh, solver, cache) - - # Iterate over all found elements. - for element in 1:n_points_curve - min_coordinate = original_nodes[:, 1, 1, 1, element_ids[element]] - max_coordinate = original_nodes[:, n_nodes, n_nodes, n_nodes, - element_ids[element]] - element_length = max_coordinate - min_coordinate - - normalized_coordinates = (curve[:, element] - min_coordinate) / - element_length[1] * 2 .- 1 - - # Interpolate to a single point in each element. - vandermonde_x = polynomial_interpolation_matrix(nodes_in, - normalized_coordinates[1]) - vandermonde_y = polynomial_interpolation_matrix(nodes_in, - normalized_coordinates[2]) - vandermonde_z = polynomial_interpolation_matrix(nodes_in, - normalized_coordinates[3]) + + return structured + end + + # Interpolate unstructured DG data to structured data (cell-centered) + # + # This function takes DG data in an unstructured, Cartesian layout and converts it to a uniformly + # distributed 2D layout. + # + # Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may + # thus be changed in future releases. + function unstructured2structured(unstructured_data, normalized_coordinates, + levels, resolution, nvisnodes_per_level) + # Extract data shape information + n_nodes_in, _, n_elements, n_variables = size(unstructured_data) + + # Get node coordinates for DG locations on reference element + nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in) + + # Calculate interpolation vandermonde matrices for each level + max_level = length(nvisnodes_per_level) - 1 + vandermonde_per_level = [] + for l in 0:max_level + n_nodes_out = nvisnodes_per_level[l + 1] + dx = 2 / n_nodes_out + nodes_out = collect(range(-1 + dx / 2, 1 - dx / 2, length = n_nodes_out)) + push!(vandermonde_per_level, + polynomial_interpolation_matrix(nodes_in, nodes_out)) + end + + # For each element, calculate index position at which to insert data in global data structure + lower_left_index = element2index(normalized_coordinates, levels, resolution, + nvisnodes_per_level) + + # Create output data structure + structured = [Matrix{Float64}(undef, resolution, resolution) for _ in 1:n_variables] + + # For each variable, interpolate element data and store to global data structure for v in 1:n_variables - for i in 1:n_nodes - for ii in 1:n_nodes - temp_data[i, ii, element, v] = (vandermonde_z * unstructured_data[i, - ii, - :, - element_ids[element], - v])[1] - end - temp_data[i, n_nodes + 1, element, v] = (vandermonde_y * temp_data[i, - 1:n_nodes, - element, - v])[1] + # Reshape data array for use in multiply_dimensionwise function + reshaped_data = reshape(unstructured_data[:, :, :, v], 1, n_nodes_in, + n_nodes_in, n_elements) + + for element_id in 1:n_elements + # Extract level for convenience + level = levels[element_id] + + # Determine target indices + n_nodes_out = nvisnodes_per_level[level + 1] + first = lower_left_index[:, element_id] + last = first .+ (n_nodes_out - 1) + + # Interpolate data + vandermonde = vandermonde_per_level[level + 1] + structured[v][first[1]:last[1], first[2]:last[2]] .= (reshape(multiply_dimensionwise(vandermonde, + reshaped_data[:, + :, + :, + element_id]), + n_nodes_out, + n_nodes_out)) end - data_on_curve[element, v] = (vandermonde_x * temp_data[:, n_nodes + 1, - element, v])[1] end - end - - return arc_length, data_on_curve, nothing -end - -# Convert 3d unstructured data from a general mesh to 1d data at given curve. -function unstructured_3d_to_1d_curve(nodes, data, curve, slice, point, nvisnodes) - # If no curve is defined, create a axis curve. - if curve === nothing - curve = axis_curve(nodes[1, :, :, :, :], nodes[2, :, :, :, :], - nodes[3, :, :, :, :], slice, point, nvisnodes) - end - - # Set up data structure. - n_points_curve = size(curve, 2) - n_variables = size(data, 1) - data_on_curve = Array{Float64}(undef, n_points_curve, n_variables) - - # Iterate over every point on the curve and determine the solutions value at given point. - for i in 1:n_points_curve - @views data_on_curve[i, :] .= get_value_at_point(curve[:, i], nodes, data) - end - - mesh_vertices_x = nothing - - return calc_arc_length(curve), data_on_curve, mesh_vertices_x -end - -# Check if the first 'amount'-many points can still form a valid tetrahedron. -function is_valid_tetrahedron(amount, coordinates; tol = 10^-4) - a = coordinates[:, 1] - b = coordinates[:, 2] - c = coordinates[:, 3] - d = coordinates[:, 4] - if amount == 2 # If two points are the same, then no tetrahedron can be formed. - return !(isapprox(a, b; atol = tol)) - elseif amount == 3 # Check if three points are on the same line. - return !on_the_same_line(a, b, c; tol = tol) - elseif amount == 4 # Check if four points form a tetrahedron. - A = hcat(coordinates[1, :], coordinates[2, :], coordinates[3, :], - SVector(1, 1, 1, 1)) - return !isapprox(det(A), 0; atol = tol) - else # With one point a tetrahedron can always be formed. - return true - end -end - -# Check if three given 3D-points are on the same line. -function on_the_same_line(a, b, c; tol = 10^-4) - # Calculate the intersection of the a-b-axis at x=0. - if b[1] == 0 - intersect_a_b = b - else - intersect_a_b = a - b .* (a[1] / b[1]) - end - # Calculate the intersection of the a-c-axis at x=0. - if c[1] == 0 - intersect_a_c = c - else - intersect_a_c = a - c .* (a[1] / c[1]) - end - return isapprox(intersect_a_b, intersect_a_c; atol = tol) -end - -# Interpolate from four corners of a tetrahedron to a single point. -function tetrahedron_interpolation(x_coordinates_in, y_coordinates_in, z_coordinates_in, - values_in, coordinate_out) - A = hcat(x_coordinates_in, y_coordinates_in, z_coordinates_in, SVector(1, 1, 1, 1)) - c = A \ values_in - return c[1] * coordinate_out[1] + c[2] * coordinate_out[2] + - c[3] * coordinate_out[3] + c[4] -end - -# Calculate the distances from every entry in node to given point. -function distances_from_single_point(nodes, point) - _, n_nodes, _, _, n_elements = size(nodes) - shifted_data = nodes .- point - distances = zeros(n_nodes, n_nodes, n_nodes, n_elements) - - # Iterate over every entry. - for element in 1:n_elements - for x in 1:n_nodes - for y in 1:n_nodes - for z in 1:n_nodes - distances[x, y, z, element] = norm(shifted_data[:, x, y, z, - element]) + + return structured + end + + # For a given normalized element coordinate, return the index of its lower left + # contribution to the global data structure + # + # Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may + # thus be changed in future releases. + function element2index3D(normalized_coordinates, levels, resolution, nvisnodes_per_level) + @assert size(normalized_coordinates, 1)==3 "only works in 3D" + + n_elements = length(levels) + + # First, determine lower left front coordinate for all cells + dx = 2 / resolution + ndim = 3 + lower_left_front_coordinate = Array{Float64}(undef, ndim, n_elements) + for element_id in 1:n_elements + nvisnodes = nvisnodes_per_level[levels[element_id] + 1] + lower_left_front_coordinate[1, element_id] = (normalized_coordinates[1, element_id] - + (nvisnodes - 1) / 2 * dx) + lower_left_front_coordinate[2, element_id] = (normalized_coordinates[2, element_id] - + (nvisnodes - 1) / 2 * dx) + lower_left_front_coordinate[3, element_id] = (normalized_coordinates[3, element_id] - + (nvisnodes - 1) / 2 * dx) + end + + # Then, convert coordinate to global index + indices = coordinate2index3D(lower_left_front_coordinate, resolution) + + return indices + end + + # Find 3D array index for a 3-tuple of normalized, cell-centered coordinates (i.e., in [-1,1]) + # + # Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may + # thus be changed in future releases. + function coordinate2index3D(coordinate, resolution::Integer) + # Calculate 1D normalized coordinates + dx = 2 / resolution + mesh_coordinates = collect(range(-1 + dx / 2, 1 - dx / 2, length = resolution)) + + # Find index + id_x = searchsortedfirst.(Ref(mesh_coordinates), coordinate[1, :], + lt = (x, y) -> x .< y .- dx / 2) + id_y = searchsortedfirst.(Ref(mesh_coordinates), coordinate[2, :], + lt = (x, y) -> x .< y .- dx / 2) + id_z = searchsortedfirst.(Ref(mesh_coordinates), coordinate[3, :], + lt = (x, y) -> x .< y .- dx / 2) + return transpose(hcat(id_x, id_y,id_z)) + end + + # For a given normalized element coordinate, return the index of its lower left + # contribution to the global data structure + # + # Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may + # thus be changed in future releases. + function element2index(normalized_coordinates, levels, resolution, nvisnodes_per_level) + @assert size(normalized_coordinates, 1)==2 "only works in 2D" + + n_elements = length(levels) + + # First, determine lower left coordinate for all cells + dx = 2 / resolution + ndim = 2 + lower_left_coordinate = Array{Float64}(undef, ndim, n_elements) + for element_id in 1:n_elements + nvisnodes = nvisnodes_per_level[levels[element_id] + 1] + lower_left_coordinate[1, element_id] = (normalized_coordinates[1, element_id] - + (nvisnodes - 1) / 2 * dx) + lower_left_coordinate[2, element_id] = (normalized_coordinates[2, element_id] - + (nvisnodes - 1) / 2 * dx) + end + + # Then, convert coordinate to global index + indices = coordinate2index(lower_left_coordinate, resolution) + + return indices + end + + # Find 2D array index for a 2-tuple of normalized, cell-centered coordinates (i.e., in [-1,1]) + # + # Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may + # thus be changed in future releases. + function coordinate2index(coordinate, resolution::Integer) + # Calculate 1D normalized coordinates + dx = 2 / resolution + mesh_coordinates = collect(range(-1 + dx / 2, 1 - dx / 2, length = resolution)) + + # Find index + id_x = searchsortedfirst.(Ref(mesh_coordinates), coordinate[1, :], + lt = (x, y) -> x .< y .- dx / 2) + id_y = searchsortedfirst.(Ref(mesh_coordinates), coordinate[2, :], + lt = (x, y) -> x .< y .- dx / 2) + return transpose(hcat(id_x, id_y)) + end + + # Calculate the vertices for each mesh cell such that it can be visualized as a closed box + # + # Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may + # thus be changed in future releases. + function calc_vertices3D(coordinates, levels, length_level_0) + ndim = size(coordinates, 1) + + # Initialize output arrays + n_elements = length(levels) + n_points_per_element = 2^ndim + 2 + x = Vector{Float64}(undef, n_points_per_element * n_elements) + y = Vector{Float64}(undef, n_points_per_element * n_elements) + z = Vector{Float64}(undef, n_points_per_element * n_elements) + + # Calculate vertices for all coordinates at once + for element_id in 1:n_elements + length = length_level_0 / 2^levels[element_id] + index = n_points_per_element * (element_id - 1) + x[index + 1] = coordinates[1, element_id] - 1 / 2 * length + x[index + 2] = coordinates[1, element_id] + 1 / 2 * length + x[index + 3] = coordinates[1, element_id] + 1 / 2 * length + x[index + 4] = coordinates[1, element_id] - 1 / 2 * length + x[index + 5] = coordinates[1, element_id] - 1 / 2 * length + x[index + 6] = coordinates[1, element_id] + 1 / 2 * length + x[index + 7] = coordinates[1, element_id] + 1 / 2 * length + x[index + 8] = coordinates[1, element_id] - 1 / 2 * length + x[index + 9] = coordinates[1, element_id] - 1 / 2 * length + x[index + 10] = NaN + + y[index + 1] = coordinates[2, element_id] - 1 / 2 * length + y[index + 2] = coordinates[2, element_id] - 1 / 2 * length + y[index + 3] = coordinates[2, element_id] + 1 / 2 * length + y[index + 4] = coordinates[2, element_id] + 1 / 2 * length + y[index + 5] = coordinates[2, element_id] - 1 / 2 * length + y[index + 6] = coordinates[2, element_id] - 1 / 2 * length + y[index + 7] = coordinates[2, element_id] + 1 / 2 * length + y[index + 8] = coordinates[2, element_id] + 1 / 2 * length + y[index + 9] = coordinates[2, element_id] - 1 / 2 * length + y[index + 10] = NaN + + z[index + 1] = coordinates[3, element_id] - 1 / 2 * length + z[index + 2] = coordinates[3, element_id] - 1 / 2 * length + z[index + 3] = coordinates[3, element_id] - 1 / 2 * length + z[index + 4] = coordinates[3, element_id] - 1 / 2 * length + z[index + 5] = coordinates[3, element_id] + 1 / 2 * length + z[index + 6] = coordinates[3, element_id] + 1 / 2 * length + z[index + 7] = coordinates[3, element_id] + 1 / 2 * length + z[index + 8] = coordinates[3, element_id] + 1 / 2 * length + z[index + 9] = coordinates[3, element_id] - 1 / 2 * length + z[index + 10] = NaN + end + + return x, y, z + end + + # Calculate the vertices to plot each grid line for StructuredMesh + # + # Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may + # thus be changed in future releases. + function calc_vertices3D(node_coordinates, mesh) + #this function uses pieces of ai-generated code + @unpack cells_per_dimension = mesh + + @assert size(node_coordinates, 1) == 3 "only works in 3D" + + linear_indices = LinearIndices(size(mesh)) + + # Initialize output arrays + n_lines = ((cells_per_dimension[2] + 1) * (cells_per_dimension[3] + 1) + # x-direction + (cells_per_dimension[1] + 1) * (cells_per_dimension[3] + 1) + # y-direction + (cells_per_dimension[1] + 1) * (cells_per_dimension[2] + 1)) # z-direction + max_length = maximum(cells_per_dimension) + n_nodes = size(node_coordinates, 2) + + x = fill(NaN, max_length * (n_nodes - 1) + 1, n_lines) + y = fill(NaN, max_length * (n_nodes - 1) + 1, n_lines) + z = fill(NaN, max_length * (n_nodes - 1) + 1, n_lines) + + line_index = 1 + + # Lines in x-direction + for cell_y in axes(mesh, 2), cell_z in axes(mesh, 3) + i = 1 + for cell_x in axes(mesh, 1) + for node in 1:(n_nodes - 1) + idx = linear_indices[cell_x, cell_y, cell_z] + x[i, line_index] = node_coordinates[1, node, 1, idx] + y[i, line_index] = node_coordinates[2, node, 1, idx] + z[i, line_index] = node_coordinates[3, node, 1, idx] + i += 1 end end + # Last point on line + x[i, line_index] = node_coordinates[1, end, 1, linear_indices[end, cell_y, cell_z]] + y[i, line_index] = node_coordinates[2, end, 1, linear_indices[end, cell_y, cell_z]] + z[i, line_index] = node_coordinates[3, end, 1, linear_indices[end, cell_y, cell_z]] + line_index += 1 end - end - return distances -end - -# Interpolate the data on given nodes to a single value at given point. -function get_value_at_point(point, nodes, data) - # Set up ata structures. - n_variables, n_x_nodes, n_y_nodes, n_z_nodes, _ = size(data) - distances = distances_from_single_point(nodes, point) - maximum_distance = maximum(distances) - - coordinates_tetrahedron = Array{Float64, 2}(undef, 3, 4) - value_tetrahedron = Array{Float64}(undef, n_variables, 4) - - index = argmin(distances) - - # If the point sits exactly on a node, no interpolation is needed. - if nodes[:, index[1], index[2], index[3], index[4]] == point - return data[1, index[1], index[2], index[3], index[4]] - end - - @views coordinates_tetrahedron[:, 1] = nodes[:, index[1], index[2], index[3], - index[4]] - @views value_tetrahedron[:, 1] = data[:, index[1], index[2], index[3], index[4]] - - # Restrict the interpolation to the closest element only. - closest_element = index[4] - @views element_distances = distances[:, :, :, closest_element] - - # Find a tetrahedron, which is given by four corners, to interpolate from. - for i in 1:4 - # Iterate until a valid tetrahedron is found. - while true - index = argmin(element_distances) - element_distances[index[1], index[2], index[3]] = maximum_distance - - @views coordinates_tetrahedron[:, i] = nodes[:, index[1], index[2], - index[3], closest_element] - @views value_tetrahedron[:, i] = data[:, index[1], index[2], index[3], - closest_element] - - # Look for another point if current tetrahedron is not valid. - if is_valid_tetrahedron(i, coordinates_tetrahedron) - break + + # Lines in y-direction + for cell_x in axes(mesh, 1), cell_z in axes(mesh, 3) + i = 1 + for cell_y in axes(mesh, 2) + for node in 1:(n_nodes - 1) + idx = linear_indices[cell_x, cell_y, cell_z] + x[i, line_index] = node_coordinates[1, 1, node, idx] + y[i, line_index] = node_coordinates[2, 1, node, idx] + z[i, line_index] = node_coordinates[3, 1, node, idx] + i += 1 + end end + # Last point on line + x[i, line_index] = node_coordinates[1, 1, end, linear_indices[cell_x, end, cell_z]] + y[i, line_index] = node_coordinates[2, 1, end, linear_indices[cell_x, end, cell_z]] + z[i, line_index] = node_coordinates[3, 1, end, linear_indices[cell_x, end, cell_z]] + line_index += 1 end - end - - # Interpolate from tetrahedron to given point. - value_at_point = Array{Float64}(undef, n_variables) - for v in 1:n_variables - value_at_point[v] = tetrahedron_interpolation(coordinates_tetrahedron[1, :], - coordinates_tetrahedron[2, :], - coordinates_tetrahedron[3, :], - value_tetrahedron[v, :], point) - end - - return value_at_point -end - -# Convert 3d unstructured data to 1d slice and interpolate them. -function unstructured_3d_to_1d(original_nodes, unstructured_data, nvisnodes, slice, - point) - if slice === :x - slice_dimension = 1 - other_dimensions = [2, 3] - elseif slice === :y - slice_dimension = 2 - other_dimensions = [1, 3] - elseif slice === :z - slice_dimension = 3 - other_dimensions = [1, 2] - else - error("illegal dimension '$slice', supported dimensions are :x, :y and :z") - end - - # Set up data structures to store new 1D data. - @views new_unstructured_data = similar(unstructured_data[1, 1, ..]) - @views temp_unstructured_data = similar(unstructured_data[1, ..]) - @views new_nodes = similar(original_nodes[1, 1, 1, ..]) - - n_nodes_in, _, _, n_elements, n_variables = size(unstructured_data) - nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in) - - # Test if point lies in the domain. - lower_limit = original_nodes[1, 1, 1, 1, 1] - upper_limit = original_nodes[1, n_nodes_in, n_nodes_in, n_nodes_in, n_elements] - - @assert length(point)>=3 "Point must be three-dimensional." - if prod(point[other_dimensions] .< lower_limit) || - prod(point[other_dimensions] .> upper_limit) - error(string("Slice axis is outside of domain. ", - " point[$other_dimensions]=$(point[other_dimensions]) must be between $lower_limit and $upper_limit")) - end - - # Count the amount of new elements. - new_id = 0 - - # Permute dimensions so that the slice dimensions are always the in correct places for later use. - if slice === :x - original_nodes = permutedims(original_nodes, [1, 3, 4, 2, 5]) - unstructured_data = permutedims(unstructured_data, [2, 3, 1, 4, 5]) - elseif slice === :y - original_nodes = permutedims(original_nodes, [1, 2, 4, 3, 5]) - unstructured_data = permutedims(unstructured_data, [1, 3, 2, 4, 5]) - end - - # Iterate over all elements to find the ones that lie on the slice axis. - for element_id in 1:n_elements - min_coordinate = original_nodes[:, 1, 1, 1, element_id] - max_coordinate = original_nodes[:, n_nodes_in, n_nodes_in, n_nodes_in, - element_id] - element_length = max_coordinate - min_coordinate - - # Test if the element is on the slice axis. If not just continue with the next element. - if !((prod(min_coordinate[other_dimensions] .<= point[other_dimensions]) && - prod(max_coordinate[other_dimensions] .> point[other_dimensions])) || - (point[other_dimensions] == upper_limit && - prod(max_coordinate[other_dimensions] .== upper_limit))) - continue - end - - new_id += 1 - - # Construct vandermonde matrix for interpolation of each 2D element to a 1D element. - normalized_intercept = (point[other_dimensions] .- - min_coordinate[other_dimensions]) / - element_length[1] * 2 .- 1 - vandermonde_i = polynomial_interpolation_matrix(nodes_in, - normalized_intercept[1]) - vandermonde_ii = polynomial_interpolation_matrix(nodes_in, - normalized_intercept[2]) - - # Interpolate to each node of new 1D element. - for v in 1:n_variables - for i in 1:n_nodes_in - for ii in 1:n_nodes_in - temp_unstructured_data[i, ii, new_id, v] = (vandermonde_ii * unstructured_data[ii, - :, - i, - element_id, - v])[1] + + # Lines in z-direction + for cell_x in axes(mesh, 1), cell_y in axes(mesh, 2) + i = 1 + for cell_z in axes(mesh, 3) + for node in 1:(n_nodes - 1) + idx = linear_indices[cell_x, cell_y, cell_z] + x[i, line_index] = node_coordinates[1, 1, node, idx] + y[i, line_index] = node_coordinates[2, 1, node, idx] + z[i, line_index] = node_coordinates[3, 1, node, idx] + i += 1 end - new_unstructured_data[i, new_id, v] = (vandermonde_i * temp_unstructured_data[i, - :, - new_id, - v])[1] end + # Last point on line + x[i, line_index] = node_coordinates[1, 1, end, linear_indices[cell_x, cell_y, end]] + y[i, line_index] = node_coordinates[2, 1, end, linear_indices[cell_x, cell_y, end]] + z[i, line_index] = node_coordinates[3, 1, end, linear_indices[cell_x, cell_y, end]] + line_index += 1 end - - new_nodes[:, new_id] = original_nodes[slice_dimension, 1, 1, :, element_id] - end - - return get_data_1d(reshape(new_nodes[:, 1:new_id], 1, n_nodes_in, new_id), - new_unstructured_data[:, 1:new_id, :], nvisnodes) -end - -# Interpolate unstructured DG data to structured data (cell-centered) -# -# This function takes DG data in an unstructured, Cartesian layout and converts it to a uniformly -# distributed 2D layout. -# -# Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may -# thus be changed in future releases. -function unstructured2structured(unstructured_data, normalized_coordinates, - levels, resolution, nvisnodes_per_level) - # Extract data shape information - n_nodes_in, _, n_elements, n_variables = size(unstructured_data) - - # Get node coordinates for DG locations on reference element - nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in) - - # Calculate interpolation vandermonde matrices for each level - max_level = length(nvisnodes_per_level) - 1 - vandermonde_per_level = [] - for l in 0:max_level - n_nodes_out = nvisnodes_per_level[l + 1] - dx = 2 / n_nodes_out - nodes_out = collect(range(-1 + dx / 2, 1 - dx / 2, length = n_nodes_out)) - push!(vandermonde_per_level, - polynomial_interpolation_matrix(nodes_in, nodes_out)) - end - - # For each element, calculate index position at which to insert data in global data structure - lower_left_index = element2index(normalized_coordinates, levels, resolution, - nvisnodes_per_level) - - # Create output data structure - structured = [Matrix{Float64}(undef, resolution, resolution) for _ in 1:n_variables] - - # For each variable, interpolate element data and store to global data structure - for v in 1:n_variables - # Reshape data array for use in multiply_dimensionwise function - reshaped_data = reshape(unstructured_data[:, :, :, v], 1, n_nodes_in, - n_nodes_in, n_elements) - + + return x, y, z + end + + # Calculate the vertices for each mesh cell such that it can be visualized as a closed box + # + # Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may + # thus be changed in future releases. + function calc_vertices(coordinates, levels, length_level_0) + ndim = size(coordinates, 1) + @assert ndim==2 "only works in 2D" + + # Initialize output arrays + n_elements = length(levels) + n_points_per_element = 2^ndim + 2 + x = Vector{Float64}(undef, n_points_per_element * n_elements) + y = Vector{Float64}(undef, n_points_per_element * n_elements) + + # Calculate vertices for all coordinates at once for element_id in 1:n_elements - # Extract level for convenience - level = levels[element_id] - - # Determine target indices - n_nodes_out = nvisnodes_per_level[level + 1] - first = lower_left_index[:, element_id] - last = first .+ (n_nodes_out - 1) - - # Interpolate data - vandermonde = vandermonde_per_level[level + 1] - structured[v][first[1]:last[1], first[2]:last[2]] .= (reshape(multiply_dimensionwise(vandermonde, - reshaped_data[:, - :, - :, - element_id]), - n_nodes_out, - n_nodes_out)) - end - end - - return structured -end - -# For a given normalized element coordinate, return the index of its lower left -# contribution to the global data structure -# -# Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may -# thus be changed in future releases. -function element2index(normalized_coordinates, levels, resolution, nvisnodes_per_level) - @assert size(normalized_coordinates, 1)==2 "only works in 2D" - - n_elements = length(levels) - - # First, determine lower left coordinate for all cells - dx = 2 / resolution - ndim = 2 - lower_left_coordinate = Array{Float64}(undef, ndim, n_elements) - for element_id in 1:n_elements - nvisnodes = nvisnodes_per_level[levels[element_id] + 1] - lower_left_coordinate[1, element_id] = (normalized_coordinates[1, element_id] - - (nvisnodes - 1) / 2 * dx) - lower_left_coordinate[2, element_id] = (normalized_coordinates[2, element_id] - - (nvisnodes - 1) / 2 * dx) - end - - # Then, convert coordinate to global index - indices = coordinate2index(lower_left_coordinate, resolution) - - return indices -end - -# Find 2D array index for a 2-tuple of normalized, cell-centered coordinates (i.e., in [-1,1]) -# -# Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may -# thus be changed in future releases. -function coordinate2index(coordinate, resolution::Integer) - # Calculate 1D normalized coordinates - dx = 2 / resolution - mesh_coordinates = collect(range(-1 + dx / 2, 1 - dx / 2, length = resolution)) - - # Find index - id_x = searchsortedfirst.(Ref(mesh_coordinates), coordinate[1, :], - lt = (x, y) -> x .< y .- dx / 2) - id_y = searchsortedfirst.(Ref(mesh_coordinates), coordinate[2, :], - lt = (x, y) -> x .< y .- dx / 2) - return transpose(hcat(id_x, id_y)) -end - -# Calculate the vertices for each mesh cell such that it can be visualized as a closed box -# -# Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may -# thus be changed in future releases. -function calc_vertices(coordinates, levels, length_level_0) - ndim = size(coordinates, 1) - @assert ndim==2 "only works in 2D" - - # Initialize output arrays - n_elements = length(levels) - n_points_per_element = 2^ndim + 2 - x = Vector{Float64}(undef, n_points_per_element * n_elements) - y = Vector{Float64}(undef, n_points_per_element * n_elements) - - # Calculate vertices for all coordinates at once - for element_id in 1:n_elements - length = length_level_0 / 2^levels[element_id] - index = n_points_per_element * (element_id - 1) - x[index + 1] = coordinates[1, element_id] - 1 / 2 * length - x[index + 2] = coordinates[1, element_id] + 1 / 2 * length - x[index + 3] = coordinates[1, element_id] + 1 / 2 * length - x[index + 4] = coordinates[1, element_id] - 1 / 2 * length - x[index + 5] = coordinates[1, element_id] - 1 / 2 * length - x[index + 6] = NaN - - y[index + 1] = coordinates[2, element_id] - 1 / 2 * length - y[index + 2] = coordinates[2, element_id] - 1 / 2 * length - y[index + 3] = coordinates[2, element_id] + 1 / 2 * length - y[index + 4] = coordinates[2, element_id] + 1 / 2 * length - y[index + 5] = coordinates[2, element_id] - 1 / 2 * length - y[index + 6] = NaN - end - - return x, y -end - -# Calculate the vertices to plot each grid line for StructuredMesh -# -# Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may -# thus be changed in future releases. -function calc_vertices(node_coordinates, mesh) - @unpack cells_per_dimension = mesh - @assert size(node_coordinates, 1)==2 "only works in 2D" - - linear_indices = LinearIndices(size(mesh)) - - # Initialize output arrays - n_lines = sum(cells_per_dimension) + 2 - max_length = maximum(cells_per_dimension) - n_nodes = size(node_coordinates, 2) - - # Create output as two matrices `x` and `y`, each holding the node locations for each of the `n_lines` grid lines - # The # of rows in the matrices must be sufficient to store the longest dimension (`max_length`), - # and for each the node locations without doubling the corner nodes (`n_nodes-1`), plus the final node (`+1`) - # Rely on Plots.jl to ignore `NaN`s (i.e., they are not plotted) to handle shorter lines - x = fill(NaN, max_length * (n_nodes - 1) + 1, n_lines) - y = fill(NaN, max_length * (n_nodes - 1) + 1, n_lines) - - line_index = 1 - # Lines in x-direction - # Bottom boundary - i = 1 - for cell_x in axes(mesh, 1) - for node in 1:(n_nodes - 1) - x[i, line_index] = node_coordinates[1, node, 1, linear_indices[cell_x, 1]] - y[i, line_index] = node_coordinates[2, node, 1, linear_indices[cell_x, 1]] - - i += 1 + length = length_level_0 / 2^levels[element_id] + index = n_points_per_element * (element_id - 1) + x[index + 1] = coordinates[1, element_id] - 1 / 2 * length + x[index + 2] = coordinates[1, element_id] + 1 / 2 * length + x[index + 3] = coordinates[1, element_id] + 1 / 2 * length + x[index + 4] = coordinates[1, element_id] - 1 / 2 * length + x[index + 5] = coordinates[1, element_id] - 1 / 2 * length + x[index + 6] = NaN + + y[index + 1] = coordinates[2, element_id] - 1 / 2 * length + y[index + 2] = coordinates[2, element_id] - 1 / 2 * length + y[index + 3] = coordinates[2, element_id] + 1 / 2 * length + y[index + 4] = coordinates[2, element_id] + 1 / 2 * length + y[index + 5] = coordinates[2, element_id] - 1 / 2 * length + y[index + 6] = NaN end - end - # Last point on bottom boundary - x[i, line_index] = node_coordinates[1, end, 1, linear_indices[end, 1]] - y[i, line_index] = node_coordinates[2, end, 1, linear_indices[end, 1]] - - # Other lines in x-direction - line_index += 1 - for cell_y in axes(mesh, 2) + + return x, y + end + + # Calculate the vertices to plot each grid line for StructuredMesh + # + # Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may + # thus be changed in future releases. + function calc_vertices(node_coordinates, mesh) + @unpack cells_per_dimension = mesh + @assert size(node_coordinates, 1)==2 "only works in 2D" + + linear_indices = LinearIndices(size(mesh)) + + # Initialize output arrays + n_lines = sum(cells_per_dimension) + 2 + max_length = maximum(cells_per_dimension) + n_nodes = size(node_coordinates, 2) + + # Create output as two matrices `x` and `y`, each holding the node locations for each of the `n_lines` grid lines + # The # of rows in the matrices must be sufficient to store the longest dimension (`max_length`), + # and for each the node locations without doubling the corner nodes (`n_nodes-1`), plus the final node (`+1`) + # Rely on Plots.jl to ignore `NaN`s (i.e., they are not plotted) to handle shorter lines + x = fill(NaN, max_length * (n_nodes - 1) + 1, n_lines) + y = fill(NaN, max_length * (n_nodes - 1) + 1, n_lines) + + line_index = 1 + # Lines in x-direction + # Bottom boundary i = 1 for cell_x in axes(mesh, 1) for node in 1:(n_nodes - 1) - x[i, line_index] = node_coordinates[1, node, end, - linear_indices[cell_x, cell_y]] - y[i, line_index] = node_coordinates[2, node, end, - linear_indices[cell_x, cell_y]] - + x[i, line_index] = node_coordinates[1, node, 1, linear_indices[cell_x, 1]] + y[i, line_index] = node_coordinates[2, node, 1, linear_indices[cell_x, 1]] + i += 1 end end - # Last point on line - x[i, line_index] = node_coordinates[1, end, end, linear_indices[end, cell_y]] - y[i, line_index] = node_coordinates[2, end, end, linear_indices[end, cell_y]] - + # Last point on bottom boundary + x[i, line_index] = node_coordinates[1, end, 1, linear_indices[end, 1]] + y[i, line_index] = node_coordinates[2, end, 1, linear_indices[end, 1]] + + # Other lines in x-direction line_index += 1 - end - - # Lines in y-direction - # Left boundary - i = 1 - for cell_y in axes(mesh, 2) - for node in 1:(n_nodes - 1) - x[i, line_index] = node_coordinates[1, 1, node, linear_indices[1, cell_y]] - y[i, line_index] = node_coordinates[2, 1, node, linear_indices[1, cell_y]] - - i += 1 + for cell_y in axes(mesh, 2) + i = 1 + for cell_x in axes(mesh, 1) + for node in 1:(n_nodes - 1) + x[i, line_index] = node_coordinates[1, node, end, + linear_indices[cell_x, cell_y]] + y[i, line_index] = node_coordinates[2, node, end, + linear_indices[cell_x, cell_y]] + + i += 1 + end + end + # Last point on line + x[i, line_index] = node_coordinates[1, end, end, linear_indices[end, cell_y]] + y[i, line_index] = node_coordinates[2, end, end, linear_indices[end, cell_y]] + + line_index += 1 end - end - # Last point on left boundary - x[i, line_index] = node_coordinates[1, 1, end, linear_indices[1, end]] - y[i, line_index] = node_coordinates[2, 1, end, linear_indices[1, end]] - - # Other lines in y-direction - line_index += 1 - for cell_x in axes(mesh, 1) + + # Lines in y-direction + # Left boundary i = 1 for cell_y in axes(mesh, 2) for node in 1:(n_nodes - 1) - x[i, line_index] = node_coordinates[1, end, node, - linear_indices[cell_x, cell_y]] - y[i, line_index] = node_coordinates[2, end, node, - linear_indices[cell_x, cell_y]] - + x[i, line_index] = node_coordinates[1, 1, node, linear_indices[1, cell_y]] + y[i, line_index] = node_coordinates[2, 1, node, linear_indices[1, cell_y]] + i += 1 end end - # Last point on line - x[i, line_index] = node_coordinates[1, end, end, linear_indices[cell_x, end]] - y[i, line_index] = node_coordinates[2, end, end, linear_indices[cell_x, end]] - + # Last point on left boundary + x[i, line_index] = node_coordinates[1, 1, end, linear_indices[1, end]] + y[i, line_index] = node_coordinates[2, 1, end, linear_indices[1, end]] + + # Other lines in y-direction line_index += 1 + for cell_x in axes(mesh, 1) + i = 1 + for cell_y in axes(mesh, 2) + for node in 1:(n_nodes - 1) + x[i, line_index] = node_coordinates[1, end, node, + linear_indices[cell_x, cell_y]] + y[i, line_index] = node_coordinates[2, end, node, + linear_indices[cell_x, cell_y]] + + i += 1 + end + end + # Last point on line + x[i, line_index] = node_coordinates[1, end, end, linear_indices[cell_x, end]] + y[i, line_index] = node_coordinates[2, end, end, linear_indices[cell_x, end]] + + line_index += 1 + end + + return x, y + end + + # Convert `slice` to orientations (1 -> `x`, 2 -> `y`, 3 -> `z`) for the two axes in a 2D plot + function _get_orientations(mesh, slice) + if ndims(mesh) == 2 || (ndims(mesh) == 3 && slice === :xy) + orientation_x = 1 + orientation_y = 2 + elseif ndims(mesh) == 3 && slice === :xz + orientation_x = 1 + orientation_y = 3 + elseif ndims(mesh) == 3 && slice === :yz + orientation_x = 2 + orientation_y = 3 + else + orientation_x = 0 + orientation_y = 0 + end + return orientation_x, orientation_y + end + + # Convert `orientation` into a guide label (see also `_get_orientations`) + function _get_guide(orientation::Integer) + if orientation == 1 + return "\$x\$" + elseif orientation == 2 + return "\$y\$" + elseif orientation == 3 + return "\$z\$" + else + return "" + end end - - return x, y -end - -# Convert `slice` to orientations (1 -> `x`, 2 -> `y`, 3 -> `z`) for the two axes in a 2D plot -function _get_orientations(mesh, slice) - if ndims(mesh) == 2 || (ndims(mesh) == 3 && slice === :xy) - orientation_x = 1 - orientation_y = 2 - elseif ndims(mesh) == 3 && slice === :xz - orientation_x = 1 - orientation_y = 3 - elseif ndims(mesh) == 3 && slice === :yz - orientation_x = 2 - orientation_y = 3 - else - orientation_x = 0 - orientation_y = 0 - end - return orientation_x, orientation_y -end - -# Convert `orientation` into a guide label (see also `_get_orientations`) -function _get_guide(orientation::Integer) - if orientation == 1 - return "\$x\$" - elseif orientation == 2 - return "\$y\$" - elseif orientation == 3 - return "\$z\$" - else - return "" - end -end - -# plotting_interpolation_matrix(dg; kwargs...) -# -# Interpolation matrix which maps discretization nodes to a set of plotting nodes. -# Defaults to the identity matrix of size `length(solver.basis.nodes)`, and interpolates -# to equispaced nodes for DGSEM (set by kwarg `nvisnodes` in the plotting function). -# -# Example: -# ```julia -# A = plotting_interpolation_matrix(dg) -# A * dg.basis.nodes # => vector of nodes at which to plot the solution -# ``` -# -# Note: we cannot use UniformScaling to define the interpolation matrix since we use it with `kron` -# to define a multi-dimensional interpolation matrix later. -plotting_interpolation_matrix(dg; kwargs...) = I(length(dg.basis.nodes)) - -function face_plotting_interpolation_matrix(dg::DGSEM; - nvisnodes = 2 * length(dg.basis.nodes)) - return polynomial_interpolation_matrix(dg.basis.nodes, LinRange(-1, 1, nvisnodes)) -end - -function plotting_interpolation_matrix(dg::DGSEM; - nvisnodes = 2 * length(dg.basis.nodes)) - Vp1D = polynomial_interpolation_matrix(dg.basis.nodes, LinRange(-1, 1, nvisnodes)) - # For quadrilateral elements, interpolation to plotting nodes involves applying a 1D interpolation - # operator to each line of nodes. This is equivalent to multiplying the vector containing all node - # node coordinates on an element by a Kronecker product of the 1D interpolation operator (e.g., a - # multi-dimensional interpolation operator). - return kron(Vp1D, Vp1D) -end - -function reference_node_coordinates_2d(dg::DGSEM) - @unpack nodes = dg.basis - r = vec([nodes[i] for i in eachnode(dg), j in eachnode(dg)]) - s = vec([nodes[j] for i in eachnode(dg), j in eachnode(dg)]) - return r, s -end - -# Find element and triangle ids containing coordinates given as a matrix [ndims, npoints] -function get_ids_by_coordinates!(ids, coordinates, pd) - if length(ids) != 2 * size(coordinates, 2) - throw(DimensionMismatch("storage length for element ids does not match the number of coordinates")) - end - - n_coordinates = size(coordinates, 2) - - for index in 1:n_coordinates - ids[index, :] .= find_element(coordinates[:, index], pd) - end - - return ids -end - -# Find the ids of elements and triangles containing given coordinates by using the triangulation in 'pd'. -function get_ids_by_coordinates(coordinates, pd) - ids = Matrix(undef, size(coordinates, 2), 2) - get_ids_by_coordinates!(ids, coordinates, pd) - return ids -end - -# Check if given 'point' is inside the triangle with corners corresponding to the coordinates of x and y. -function is_in_triangle(point, x, y) - a = SVector(x[1], y[1]) - b = SVector(x[2], y[2]) - c = SVector(x[3], y[3]) - return is_on_same_side(point, a, b, c) && is_on_same_side(point, b, c, a) && - is_on_same_side(point, c, a, b) -end - -# Create an axis through x and y to then check if 'point' is on the same side of the axis as z. -function is_on_same_side(point, x, y, z) - if (y[1] - x[1]) == 0 - return (point[1] - x[1]) * (z[1] - x[1]) >= 0 - else - a = (y[2] - x[2]) / (y[1] - x[1]) - b = x[2] - a * x[1] - return (z[2] - a * z[1] - b) * (point[2] - a * point[1] - b) >= 0 + + # plotting_interpolation_matrix(dg; kwargs...) + # + # Interpolation matrix which maps discretization nodes to a set of plotting nodes. + # Defaults to the identity matrix of size `length(solver.basis.nodes)`, and interpolates + # to equispaced nodes for DGSEM (set by kwarg `nvisnodes` in the plotting function). + # + # Example: + # ```julia + # A = plotting_interpolation_matrix(dg) + # A * dg.basis.nodes # => vector of nodes at which to plot the solution + # ``` + # + # Note: we cannot use UniformScaling to define the interpolation matrix since we use it with `kron` + # to define a multi-dimensional interpolation matrix later. + plotting_interpolation_matrix(dg; kwargs...) = I(length(dg.basis.nodes)) + + function face_plotting_interpolation_matrix(dg::DGSEM; + nvisnodes = 2 * length(dg.basis.nodes)) + return polynomial_interpolation_matrix(dg.basis.nodes, LinRange(-1, 1, nvisnodes)) + end + + function plotting_interpolation_matrix(dg::DGSEM; + nvisnodes = 2 * length(dg.basis.nodes)) + Vp1D = polynomial_interpolation_matrix(dg.basis.nodes, LinRange(-1, 1, nvisnodes)) + # For quadrilateral elements, interpolation to plotting nodes involves applying a 1D interpolation + # operator to each line of nodes. This is equivalent to multiplying the vector containing all node + # node coordinates on an element by a Kronecker product of the 1D interpolation operator (e.g., a + # multi-dimensional interpolation operator). + return kron(Vp1D, Vp1D) + end + + function reference_node_coordinates_2d(dg::DGSEM) + @unpack nodes = dg.basis + r = vec([nodes[i] for i in eachnode(dg), j in eachnode(dg)]) + s = vec([nodes[j] for i in eachnode(dg), j in eachnode(dg)]) + return r, s + end + + # Find element and triangle ids containing coordinates given as a matrix [ndims, npoints] + function get_ids_by_coordinates!(ids, coordinates, pd) + if length(ids) != 2 * size(coordinates, 2) + throw(DimensionMismatch("storage length for element ids does not match the number of coordinates")) + end + + n_coordinates = size(coordinates, 2) + + for index in 1:n_coordinates + ids[index, :] .= find_element(coordinates[:, index], pd) + end + + return ids + end + + # Find the ids of elements and triangles containing given coordinates by using the triangulation in 'pd'. + function get_ids_by_coordinates(coordinates, pd) + ids = Matrix(undef, size(coordinates, 2), 2) + get_ids_by_coordinates!(ids, coordinates, pd) + return ids + end + + # Check if given 'point' is inside the triangle with corners corresponding to the coordinates of x and y. + function is_in_triangle(point, x, y) + a = SVector(x[1], y[1]) + b = SVector(x[2], y[2]) + c = SVector(x[3], y[3]) + return is_on_same_side(point, a, b, c) && is_on_same_side(point, b, c, a) && + is_on_same_side(point, c, a, b) + end + + # Create an axis through x and y to then check if 'point' is on the same side of the axis as z. + function is_on_same_side(point, x, y, z) + if (y[1] - x[1]) == 0 + return (point[1] - x[1]) * (z[1] - x[1]) >= 0 + else + a = (y[2] - x[2]) / (y[1] - x[1]) + b = x[2] - a * x[1] + return (z[2] - a * z[1] - b) * (point[2] - a * point[1] - b) >= 0 + end end -end - -# For a given 'point', return the id of the element it is contained in in; if not found return 0. -function find_element(point, pd) - n_tri = size(pd.t, 1) - n_elements = size(pd.x, 2) - - # Iterate over all elements. - for element in 1:n_elements - # Iterate over all triangles in given element. - for tri in 1:n_tri - if is_in_triangle(point, pd.x[pd.t[tri, :], element], - pd.y[pd.t[tri, :], element]) - return SVector(element, tri) + + # For a given 'point', return the id of the element it is contained in in; if not found return 0. + function find_element(point, pd) + n_tri = size(pd.t, 1) + n_elements = size(pd.x, 2) + + # Iterate over all elements. + for element in 1:n_elements + # Iterate over all triangles in given element. + for tri in 1:n_tri + if is_in_triangle(point, pd.x[pd.t[tri, :], element], + pd.y[pd.t[tri, :], element]) + return SVector(element, tri) + end end end end -end - -# Interpolate from three corners of a triangle to a single point. -function triangle_interpolation(x_coordinates_in, y_coordinates_in, values_in, - coordinate_out) - A = hcat(x_coordinates_in, y_coordinates_in, SVector(1, 1, 1)) - c = A \ values_in - return c[1] * coordinate_out[1] + c[2] * coordinate_out[2] + c[3] -end - -# Create an axis. -function axis_curve(nodes_x, nodes_y, nodes_z, slice, point, n_points) - if n_points === nothing - n_points = 64 - end - dimensions = length(point) - curve = zeros(dimensions, n_points) - if slice == :x - xmin, xmax = extrema(nodes_x) - curve[1, :] .= range(xmin, xmax, length = n_points) - curve[2, :] .= point[2] - if dimensions === 3 - curve[3, :] .= point[3] - end - elseif slice == :y - ymin, ymax = extrema(nodes_y) - curve[1, :] .= point[1] - curve[2, :] .= range(ymin, ymax, length = n_points) - if dimensions === 3 - curve[3, :] .= point[3] - end - elseif slice == :z - zmin, zmax = extrema(nodes_z) - curve[1, :] .= point[1] - curve[2, :] .= point[2] - curve[3, :] .= range(zmin, zmax, length = n_points) - else - @assert false "Input for 'slice' is not supported here." + + # Interpolate from three corners of a triangle to a single point. + function triangle_interpolation(x_coordinates_in, y_coordinates_in, values_in, + coordinate_out) + A = hcat(x_coordinates_in, y_coordinates_in, SVector(1, 1, 1)) + c = A \ values_in + return c[1] * coordinate_out[1] + c[2] * coordinate_out[2] + c[3] + end + + # Create an axis. + function axis_curve(nodes_x, nodes_y, nodes_z, slice, point, n_points) + if n_points === nothing + n_points = 64 + end + dimensions = length(point) + curve = zeros(dimensions, n_points) + if slice == :x + xmin, xmax = extrema(nodes_x) + curve[1, :] .= range(xmin, xmax, length = n_points) + curve[2, :] .= point[2] + if dimensions === 3 + curve[3, :] .= point[3] + end + elseif slice == :y + ymin, ymax = extrema(nodes_y) + curve[1, :] .= point[1] + curve[2, :] .= range(ymin, ymax, length = n_points) + if dimensions === 3 + curve[3, :] .= point[3] + end + elseif slice == :z + zmin, zmax = extrema(nodes_z) + curve[1, :] .= point[1] + curve[2, :] .= point[2] + curve[3, :] .= range(zmin, zmax, length = n_points) + else + @assert false "Input for 'slice' is not supported here." + end + + return curve end - - return curve -end -end # @muladd + end # @muladd \ No newline at end of file From 23b5f2f0c4a27e7e5021aab8040b56a0cbfc1f24 Mon Sep 17 00:00:00 2001 From: s6nistam Date: Fri, 13 Dec 2024 07:41:18 +0100 Subject: [PATCH 4/6] visualization Callback using Makie removing GLMakie and Plots from project toml --- Project.toml | 2 - ...ixir_euler_blast_wave_amr_visualization.jl | 109 ++++++++++++++++++ src/callbacks_step/visualization.jl | 14 --- 3 files changed, 109 insertions(+), 16 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_euler_blast_wave_amr_visualization.jl diff --git a/Project.toml b/Project.toml index ba3f3117f47..8a3a4772613 100644 --- a/Project.toml +++ b/Project.toml @@ -15,7 +15,6 @@ Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -26,7 +25,6 @@ MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Preferences = "21216c6a-2e73-6563-6e65-726566657250" diff --git a/examples/tree_2d_dgsem/elixir_euler_blast_wave_amr_visualization.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_amr_visualization.jl new file mode 100644 index 00000000000..b4863e4c382 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_blast_wave_amr_visualization.jl @@ -0,0 +1,109 @@ + +using OrdinaryDiffEq +using Trixi +using GLMakie + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +""" + initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) + +A medium blast wave taken from +- Sebastian Hennemann, Gregor J. Gassner (2020) + A provably entropy stable subcell shock capturing approach for high order split form DG + [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) +""" +function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) + # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" + # Set up polar coordinates + inicenter = SVector(0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + phi = atan(y_norm, x_norm) + sin_phi, cos_phi = sincos(phi) + + # Calculate primitive variables + rho = r > 0.5 ? 1.0 : 1.1691 + v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi + v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi + p = r > 0.5 ? 1.0E-3 : 1.245 + + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_blast_wave + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +basis = LobattoLegendreBasis(3) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-2.0, -2.0) +coordinates_max = (2.0, 2.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 6, + n_cells_max = 10_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 12.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +amr_indicator = IndicatorHennemannGassner(semi, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) + + +visualization = VisualizationCallback(interval = 20, clims = (0, 1), plot_creator = Trixi.show_plot_makie) + +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 4, + max_level = 6, max_threshold = 0.01) +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +stepsize_callback = StepsizeCallback(cfl = 0.9) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + visualization, + amr_callback, stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/callbacks_step/visualization.jl b/src/callbacks_step/visualization.jl index 7d34e63f16e..c95266c860d 100644 --- a/src/callbacks_step/visualization.jl +++ b/src/callbacks_step/visualization.jl @@ -102,20 +102,6 @@ function VisualizationCallback(; interval = 0, plot_data_creator, plot_creator, nothing, Dict{Symbol, Any}(plot_arguments)) - # Warn users if they create a visualization callback without having loaded the Plots package - # - # Note: This warning is added for convenience, as Plots is the only "officially" supported - # visualization package right now. However, in general nothing prevents anyone from using - # other packages such as Makie, Gadfly etc., given that appropriate `plot_creator`s are - # passed. This is also the reason why the visualization callback is not included via - # Requires.jl only when Plots is present. - # In the future, we should update/remove this warning if other plotting packages are - # starting to be used. - if !(:Plots in names(@__MODULE__, all = true) || :Makie in names(@__MODULE__, all = true)) - println(names(@__MODULE__, all = true)) - @warn "Package `Plots` or `Makie` required by `VisualizationCallback` to visualize results" - end - DiscreteCallback(visualization_callback, visualization_callback, # the first one is the condition, the second the affect! save_positions = (false, false), initialize = initialize!) From 5866ae71e81a31050f7295b7ddd26ff75d22c4f8 Mon Sep 17 00:00:00 2001 From: s6nistam Date: Fri, 13 Dec 2024 08:07:33 +0100 Subject: [PATCH 5/6] Visualization Callback Makie added Labels --- src/callbacks_step/visualization.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/callbacks_step/visualization.jl b/src/callbacks_step/visualization.jl index c95266c860d..8ebe69e0ca7 100644 --- a/src/callbacks_step/visualization.jl +++ b/src/callbacks_step/visualization.jl @@ -297,11 +297,11 @@ function show_plot_makie(ndims, visualization_callback, plot_data, variable_name figure = visualization_callback.figure if ndims == 2 for v in 1:size(variable_names)[1] - GLMakie.heatmap(figure[makieLayoutHelper(v)...], plot_data.x, plot_data.y, plot_data.data[v]) + GLMakie.heatmap(figure[makieLayoutHelper(v)...], plot_data.x, plot_data.y, plot_data.data[v], axis=(; title = variable_names[v])) end else for v in 1:size(variable_names)[1] - GLMakie.volume(figure[makieLayoutHelper(v)...], plot_data.data[v]) + GLMakie.volume(figure[makieLayoutHelper(v)...], plot_data.data[v], axis=(; title = variable_names[v])) end end GLMakie.display(figure) From a1ce5ef82ceca09bd31b2c5e9c2862e12ba903ab Mon Sep 17 00:00:00 2001 From: s6nistam Date: Fri, 13 Dec 2024 09:52:59 +0100 Subject: [PATCH 6/6] Visualization Callback with makie refactored to use axis to fix duplicate issue and to make labels in 3D work --- .../elixir_euler_amr_visualization.jl | 87 +++++++++++++++++++ src/callbacks_step/visualization.jl | 28 +++--- 2 files changed, 99 insertions(+), 16 deletions(-) create mode 100644 examples/tree_3d_dgsem/elixir_euler_amr_visualization.jl diff --git a/examples/tree_3d_dgsem/elixir_euler_amr_visualization.jl b/examples/tree_3d_dgsem/elixir_euler_amr_visualization.jl new file mode 100644 index 00000000000..643ea465db7 --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_euler_amr_visualization.jl @@ -0,0 +1,87 @@ + +using OrdinaryDiffEq +using Trixi +using Plots +using GLMakie + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +""" + initial_condition_density_pulse(x, t, equations::CompressibleEulerEquations3D) + +A Gaussian pulse in the density with constant velocity and pressure; reduces the +compressible Euler equations to the linear advection equations. +""" +function initial_condition_density_pulse(x, t, equations::CompressibleEulerEquations3D) + rho = 1 + exp(-(x[1]^2 + x[2]^2 + x[3]^2)) / 2 + v1 = 1 + v2 = 1 + v3 = 1 + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_v3 = rho * v3 + p = 1 + rho_e = p / (equations.gamma - 1) + 1 / 2 * rho * (v1^2 + v2^2 + v3^2) + return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e) +end +initial_condition = initial_condition_density_pulse +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + +coordinates_min = (-5.0, -5.0, -5.0) +coordinates_max = (5.0, 5.0, 5.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 10_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (entropy,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_restart = SaveRestartCallback(interval = 100, + save_final_restart = true) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +visualization = VisualizationCallback(interval = 20, clims = (0, 1), plot_data_creator = Trixi.PlotData3D, plot_creator = Trixi.show_plot_makie) + +amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first), + base_level = 1, + med_level = 2, med_threshold = 1.05, + max_level = 3, max_threshold = 1.3) +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +stepsize_callback = StepsizeCallback(cfl = 0.9) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_restart, save_solution, visualization, + amr_callback, stepsize_callback); + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/callbacks_step/visualization.jl b/src/callbacks_step/visualization.jl index 8ebe69e0ca7..e03128ef726 100644 --- a/src/callbacks_step/visualization.jl +++ b/src/callbacks_step/visualization.jl @@ -14,6 +14,7 @@ mutable struct VisualizationCallback{SolutionVariables, VariableNames, PlotDataC plot_data_creator::PlotDataCreator plot_creator::PlotCreator figure + axis plot_arguments::Dict{Symbol, Any} end @@ -99,7 +100,7 @@ function VisualizationCallback(; interval = 0, visualization_callback = VisualizationCallback(interval, solution_variables, variable_names, show_mesh, - plot_data_creator, plot_creator, nothing, + plot_data_creator, plot_creator, nothing, [], Dict{Symbol, Any}(plot_arguments)) DiscreteCallback(visualization_callback, visualization_callback, # the first one is the condition, the second the affect! @@ -137,7 +138,7 @@ end function (visualization_callback::VisualizationCallback)(integrator) u_ode = integrator.u semi = integrator.p - @unpack plot_arguments, solution_variables, variable_names, show_mesh, plot_data_creator, plot_creator, figure = visualization_callback + @unpack plot_arguments, solution_variables, variable_names, show_mesh, plot_data_creator, plot_creator, figure, axis = visualization_callback mesh, equations, solver, cache = mesh_equations_solver_cache(integrator.p) n = Trixi.ndims(mesh) @@ -152,7 +153,7 @@ function (visualization_callback::VisualizationCallback)(integrator) # Create plot plot_creator(n, visualization_callback, plot_data, variable_names; show_mesh = show_mesh, plot_arguments = plot_arguments, - time = integrator.t, timestep = integrator.stats.naccept, figure = figure) + time = integrator.t, timestep = integrator.stats.naccept, figure = figure, axis = axis) # avoid re-evaluating possible FSAL stages u_modified!(integrator, false) @@ -178,7 +179,7 @@ See also: [`VisualizationCallback`](@ref), [`save_plot`](@ref) """ function show_plot(ndims, visualization_callback, plot_data, variable_names; show_mesh = true, plot_arguments = Dict{Symbol, Any}(), - time = nothing, timestep = nothing, figure = nothing) + time = nothing, timestep = nothing, figure = nothing, axis = nothing) # Gather subplots plots = [] for v in variable_names @@ -289,30 +290,25 @@ See also: [`VisualizationCallback`](@ref), [`save_plot`](@ref) """ function show_plot_makie(ndims, visualization_callback, plot_data, variable_names; show_mesh = true, plot_arguments = Dict{Symbol, Any}(), - time = nothing, timestep = nothing, figure = nothing) - + time = nothing, timestep = nothing, figure = nothing, axis = []) if figure === nothing @warn "Creating new figure" visualization_callback.figure = GLMakie.Figure() figure = visualization_callback.figure - if ndims == 2 - for v in 1:size(variable_names)[1] - GLMakie.heatmap(figure[makieLayoutHelper(v)...], plot_data.x, plot_data.y, plot_data.data[v], axis=(; title = variable_names[v])) - end - else - for v in 1:size(variable_names)[1] - GLMakie.volume(figure[makieLayoutHelper(v)...], plot_data.data[v], axis=(; title = variable_names[v])) - end + for v in 1:size(variable_names)[1] + push!(axis, (ndims == 2) ? GLMakie.Axis(figure[makieLayoutHelper(v)...], title = variable_names[v]) : + GLMakie.Axis3(figure[makieLayoutHelper(v)...], aspect=:equal, title = variable_names[v])) end + visualization_callback.axis = axis GLMakie.display(figure) else if ndims == 2 for v in 1:size(variable_names)[1] - GLMakie.heatmap!(figure[makieLayoutHelper(v)...], plot_data.x, plot_data.y, plot_data.data[v]) + GLMakie.heatmap!(axis[v], plot_data.x, plot_data.y, plot_data.data[v]) end else for v in 1:size(variable_names)[1] - GLMakie.volume!(figure[makieLayoutHelper(v)...], plot_data.data[v]) + GLMakie.volume!(axis[v], plot_data.data[v]) end end end