Clean up VTK utilities and tests

[jlchan]

Improve type stability and simplify code

[codecov]

Codecov Report

Merging #87 (abbdb23) into main (c49b289) will increase coverage by 0.03%.
The diff coverage is 93.75%.

❗ Current head abbdb23 differs from pull request most recent head ad918b2. Consider uploading reports for the commit ad918b2 to get more accurate results

@@            Coverage Diff             @@
##             main      #87      +/-   ##
==========================================
+ Coverage   96.57%   96.60%   +0.03%     
==========================================
  Files          22       22              
  Lines        2424     2416       -8     
==========================================
- Hits         2341     2334       -7     
+ Misses         83       82       -1     

Impacted Files	Coverage Δ
src/StartUpDG.jl	100.00% <ø> (ø)
src/mesh/mesh_visualization.jl	94.33% <83.33%> (+1.11%)	⬆️
src/mesh/vtk_helper.jl	100.00% <100.00%> (ø)

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[jlchan]

MeshData_to_vtk is still type unstable, but I think that's due to JuliaVTK/WriteVTK.jl#116.

[Davknapp]

For triangles we need a return value to later compute the inner nodes for prisms. There we subtract 3 from the prism-order and call this function with it.

[jlchan]

Hm, previously for order < 0 we returned nothing

    if order < 0
        return nothing
    end

How is the return value of nothing used? It would be nice if we could avoid this, since Julia can't infer at compile time what return type to expect (which causes type instability).

[Davknapp]

I probably can fix it in the prism-branch. We can check if the order stays positiv. If it is negative we don't need to append the coordiantes of the prism-nodes by the triangular nodes.
Therefore, this assertion can stay.

[Davknapp]

Thanks for clean-up, it helps to learn all the tricks in Julia! We can merge, as soon as the triangle-order-problem is solved

[jipolanco]

MeshData_to_vtk is still type unstable, but I think that's due to jipolanco/WriteVTK.jl#116.

I haven't looked at this package's code in detail, but I think it's rather because the type of coords in MeshData_to_vtk is not inferred, which leads to vtkfile not being inferred either.

Minimal example:

using StartUpDG

# polynomial degree and mesh size
N = 3
K1D = 8

# init ref element and mesh
rd = RefElemData(Tri(), N)
VXY, EToV = uniform_mesh(Tri(), K1D)
md = MeshData(VXY, EToV, rd)

MeshData_to_vtk(md, rd, nothing, nothing, "testfile", false, true)

julia> @code_warntype MeshData_to_vtk(md, rd, nothing, nothing, "testfile", false, true)
MethodInstance for StartUpDG.MeshData_to_vtk(::MeshData{2, StartUpDG.VertexMappedMesh{Tri, Tuple{Vector{Float64}, Vector{Float64}}, Matrix{Int64}}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Int64}, Matrix{Float64}, StaticArraysCore.SMatrix{2, 2, Matrix{Float64}, 4}, Matrix{Float64}, Matrix{Int64}, Matrix{Int64}, Vector{Int64}}, ::RefElemData{2, Tri, Polynomial, Tuple{Vector{Int64}, Vector{Int64}, Vector{Int64}}, Tuple{Vector{Float64}, Vector{Float64}}, Tuple{Vector{Float64}, Vector{Float64}}, Tuple{Vector{Float64}, Vector{Float64}}, Tuple{Vector{Float64}, Vector{Float64}}, Tuple{Vector{Float64}, Vector{Float64}}, Vector{Int64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Tuple{Matrix{Float64}, Matrix{Float64}}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Vector{Float64}, Vector{Float64}}, ::Nothing, ::Nothing, ::String, ::Bool, ::Bool)
  from MeshData_to_vtk(md::MeshData, rd::RefElemData, data, dataname, filename, write_data, equi_dist_nodes) @ StartUpDG ~/.julia/dev/StartUpDG/src/mesh/mesh_visualization.jl:82
Arguments
  #self#::Core.Const(StartUpDG.MeshData_to_vtk)
  md::MeshData{2, StartUpDG.VertexMappedMesh{Tri, Tuple{Vector{Float64}, Vector{Float64}}, Matrix{Int64}}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Int64}, Matrix{Float64}, StaticArraysCore.SMatrix{2, 2, Matrix{Float64}, 4}, Matrix{Float64}, Matrix{Int64}, Matrix{Int64}, Vector{Int64}}
  rd::RefElemData{2, Tri, Polynomial, Tuple{Vector{Int64}, Vector{Int64}, Vector{Int64}}, Tuple{Vector{Float64}, Vector{Float64}}, Tuple{Vector{Float64}, Vector{Float64}}, Tuple{Vector{Float64}, Vector{Float64}}, Tuple{Vector{Float64}, Vector{Float64}}, Tuple{Vector{Float64}, Vector{Float64}}, Vector{Int64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Tuple{Matrix{Float64}, Matrix{Float64}}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Vector{Float64}, Vector{Float64}}
  data::Core.Const(nothing)
  dataname::Core.Const(nothing)
  filename::String
  write_data::Bool
  equi_dist_nodes::Bool
Locals
  @_9::Union{}
  #346::StartUpDG.var"#346#348"
  #345::StartUpDG.var"#345#347"{VTKCellType, Int64, Vector{Int64}}
  vtkfile::Any
  coords::Tuple{Any, Any}
  interpolate::Core.Box
  cells::Vector{MeshCell{VTKCellType, Vector{Int64}}}
  vtk_cell_type::VTKCellType
  num_lagrange_points::Int64
  perm::Vector{Int64}
  i::Union{}
Body::Any
1 ─       Core.NewvarNode(:(@_9))
│         Core.NewvarNode(:(#346))
│         Core.NewvarNode(:(vtkfile))
│         Core.NewvarNode(:(coords))
│         (interpolate = Core.Box())
│         (perm = StartUpDG.SUD_to_vtk_order(rd))
│         (num_lagrange_points = StartUpDG.length(perm))
│   %8  = Base.getproperty(rd, :element_type)::Core.Const(Tri())
│         (vtk_cell_type = StartUpDG.type_to_vtk(%8))
│   %10 = StartUpDG.:(var"#345#347")::Core.Const(StartUpDG.var"#345#347")
│   %11 = Core.typeof(vtk_cell_type::Core.Const(VTKCellType("VTK_LAGRANGE_TRIANGLE", 0x45, -1)))::Core.Const(VTKCellType)
│   %12 = Core.typeof(num_lagrange_points)::Core.Const(Int64)
│   %13 = Core.typeof(perm)::Core.Const(Vector{Int64})
│   %14 = Core.apply_type(%10, %11, %12, %13)::Core.Const(StartUpDG.var"#345#347"{VTKCellType, Int64, Vector{Int64}})
│   %15 = vtk_cell_type::Core.Const(VTKCellType("VTK_LAGRANGE_TRIANGLE", 0x45, -1))
│   %16 = num_lagrange_points::Int64
│         (#345 = %new(%14, %15, %16, perm))
│   %18 = #345::Core.PartialStruct(StartUpDG.var"#345#347"{VTKCellType, Int64, Vector{Int64}}, Any[Core.Const(VTKCellType("VTK_LAGRANGE_TRIANGLE", 0x45, -1)), Int64, Vector{Int64}])
│   %19 = Base.getproperty(md, :num_elements)::Int64
│   %20 = (1:%19)::Core.PartialStruct(UnitRange{Int64}, Any[Core.Const(1), Int64])
│   %21 = Base.Generator(%18, %20)::Core.PartialStruct(Base.Generator{UnitRange{Int64}, StartUpDG.var"#345#347"{VTKCellType, Int64, Vector{Int64}}}, Any[Core.PartialStruct(StartUpDG.var"#345#347"{VTKCellType, Int64, Vector{Int64}}, Any[Core.Const(VTKCellType("VTK_LAGRANGE_TRIANGLE", 0x45, -1)), Int64, Vector{Int64}]), Core.PartialStruct(UnitRange{Int64}, Any[Core.Const(1), Int64])])
│         (cells = Base.collect(%21))
│   %23 = (equi_dist_nodes == true)::Bool
└──       goto #3 if not %23
2 ─ %25 = StartUpDG.vandermonde::Core.Const(NodesAndModes.vandermonde)
│   %26 = Base.getproperty(rd, :element_type)::Core.Const(Tri())
│   %27 = Base.getproperty(rd, :N)::Int64
│   %28 = Core.tuple(%26, %27)::Tuple{Tri, Int64}
│   %29 = Base.getproperty(rd, :element_type)::Core.Const(Tri())
│   %30 = Base.getproperty(rd, :N)::Int64
│   %31 = StartUpDG.equi_nodes(%29, %30)::Tuple{Vector{Float64}, Vector{Float64}}
│   %32 = Core._apply_iterate(Base.iterate, %25, %28, %31)::Matrix{Float64}
│   %33 = Base.getproperty(rd, :VDM)::Matrix{Float64}
│   %34 = (%32 / %33)::Matrix{Float64}
│         Core.setfield!(interpolate, :contents, %34)
└──       goto #4
3 ─ %37 = StartUpDG.I(num_lagrange_points)::LinearAlgebra.Diagonal{Bool, Vector{Bool}}
└──       Core.setfield!(interpolate, :contents, %37)
4 ┄       (#346 = %new(StartUpDG.:(var"#346#348"), interpolate))
│   %40 = #346::StartUpDG.var"#346#348"
│   %41 = Base.getproperty(md, :xyz)::Tuple{Matrix{Float64}, Matrix{Float64}}
│         (coords = StartUpDG.map(%40, %41))
│   %43 = StartUpDG.vtk_grid::Core.Const(WriteVTK.vtk_grid)
│   %44 = Core.tuple(filename)::Tuple{String}
│   %45 = coords::Tuple{Any, Any}
│   %46 = Core.tuple(cells)::Tuple{Vector{MeshCell{VTKCellType, Vector{Int64}}}}
│         (vtkfile = Core._apply_iterate(Base.iterate, %43, %44, %45, %46))
└──       goto #6 if not write_data
5 ─       StartUpDG.length(dataname)
│         Core.Const(:(1:%49))
│         Core.Const(:(@_9 = Base.iterate(%50)))
│         Core.Const(:(@_9 === nothing))
│         Core.Const(:(Base.not_int(%52)))
│         Core.Const(:(goto %67 if not %53))
│         Core.Const(:(@_9))
│         Core.Const(:(i = Core.getfield(%55, 1)))
│         Core.Const(:(Core.getfield(%55, 2)))
│         Core.Const(:(Base.getindex(data, i)))
│         Core.Const(:(vtkfile))
│         Core.Const(:(Base.getindex(dataname, i)))
│         Core.Const(:(Base.setindex!(%59, %58, %60)))
│         Core.Const(:(@_9 = Base.iterate(%50, %57)))
│         Core.Const(:(@_9 === nothing))
│         Core.Const(:(Base.not_int(%63)))
│         Core.Const(:(goto %67 if not %64))
└──       Core.Const(:(goto %55))
6 ┄ %67 = StartUpDG.vtk_save(vtkfile)::Any
└──       return %67

[jlchan]

@jipolanco thanks for catching that! Looks like the splatting might be causing issues here.

[Davknapp]

In my opinion it is ready to merge!