Merge pull request #15 from ZenanH/main

🎉 Release new version v0.2.0

Author: ZenanH

Project.toml

@@ -1,7 +1,7 @@
 name = "MaterialPointVisualizer"
 uuid = "9ce2fbfb-c269-402f-8683-a675189e795c"
 authors = ["ZenanH <zenan.huo@outlook.com>"]
-version = "0.1.11"
+version = "0.2.0"
 
 [deps]
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
@@ -12,18 +12,20 @@ HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 WGLMakie = "276b4fcb-3e11-5398-bf8b-a0c2d153d008"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
-ColorSchemes = "3.30"
+ColorSchemes = "3.31"
 DelimitedFiles = "1"
-FastPointQuery = "0.1"
+FastPointQuery = "0.2"
 HDF5 = "0.17"
 PrecompileTools = "1.2"
+ProgressMeter = "1.10"
 StatsBase = "0.34"
-WGLMakie = "0.11"
+WGLMakie = "0.13"
 WriteVTK = "1.21"
 julia = "1.11"
 

README.md

@@ -3,7 +3,7 @@
 [![CI](https://github.com/LandslideSIM/MaterialPointVisualizer.jl/actions/workflows/ci.yml/badge.svg)](https://github.com/LandslideSIM/MaterialPointVisualizer.jl/actions/workflows/ci.yml) 
 [![Stable](https://img.shields.io/badge/docs-stable-blue.svg?logo=quicklook)](https://LandslideSIM.github.io/MaterialPointVisualizer.jl/stable)
 [![Dev](https://img.shields.io/badge/docs-dev-red.svg?logo=quicklook)](https://LandslideSIM.github.io/MaterialPointVisualizer.jl/dev)
-[![Version](https://img.shields.io/badge/version-v0.1.11-pink)]()
+[![Version](https://img.shields.io/badge/version-v0.2.0-pink)]()
 
 With this package, we can convert the MPM simulation results (HDF5 files from ***[MaterialPointSolver.jl](https://github.com/LandslideSIM/MaterialPointSolver.jl)*** ) into `.vtp` files or create ParaView-compatible animations. Additionally, it includes some post-processing functionalities.
 

src/MaterialPointVisualizer.jl

@@ -11,7 +11,8 @@
 
 module MaterialPointVisualizer
 
-using ColorSchemes, Dates, DelimitedFiles, HDF5, Logging, PrecompileTools, Printf, WGLMakie, WriteVTK
+using ColorSchemes, Dates, DelimitedFiles, HDF5, Logging, PrecompileTools, Printf,
+    ProgressMeter, #=WGLMakie,=# WriteVTK
 
 import StatsBase: sample
 import FastPointQuery: trimesh, splashsurf, np, meshio
@@ -26,10 +27,10 @@ import FastPointQuery: trimesh, splashsurf, np, meshio
         @sprintf("%02d days: %s", days, Dates.format(time, "HH:MM:SS"))
 end
 
-include(joinpath(@__DIR__, "hdf2pvd.jl"     ))
-include(joinpath(@__DIR__, "pts2vtp.jl"     ))
-include(joinpath(@__DIR__, "pts2surf.jl"    ))
-include(joinpath(@__DIR__, "plot/display.jl"))
+include(joinpath(@__DIR__, "hdf2pvd.jl"))
+include(joinpath(@__DIR__, "pts2vtp.jl"))
+include(joinpath(@__DIR__, "pts2surf.jl"))
+# include(joinpath(@__DIR__, "plot/display.jl"))
 
 quiet(f) = redirect_stdout(devnull) do
     redirect_stderr(devnull) do
@@ -39,6 +40,6 @@ quiet(f) = redirect_stdout(devnull) do
     end
 end
 
-include(joinpath(@__DIR__, "precompile.jl"))
+# include(joinpath(@__DIR__, "precompile.jl"))
 
 end

src/hdf2pvd.jl

@@ -39,7 +39,7 @@ If your HDF5 file is generated by other means, please refer to [WriteVTK.jl](htt
     h5_str = joinpath(conf.prjdst, "$(conf.prjname).h5")
     h5_file = isfile(h5_str) ? h5_str : throw(ArgumentError("Expected a h5 file path at $h5_str"))
     anim_path = joinpath(conf.prjdst, "animation"); mkpath(anim_path)
-    pvd = paraview_collection(joinpath(conf.prjdst, "$(conf.prjname).pvd"))
+    #pvd = paraview_collection(joinpath(conf.prjdst, "$(conf.prjname).pvd"))
     fid = h5open(h5_file, "r")
     required = ("time", "ξ")
     groups = Tuple{String,Float64}[]
@@ -51,34 +51,32 @@ If your HDF5 file is generated by other means, please refer to [WriteVTK.jl](htt
     end; sort!(groups, by = last)
     isempty(groups) && @warn "No valid time-step groups found in $h5_str"
 
-    p_total, p_iters = length(groups), 0
-    t1, t_start = time(), time()
-    @inbounds for (step, (gname, t)) in enumerate(groups)
-        grp = fid[gname]
-        coords = read(grp["ξ"])
-        vtk_path = joinpath(anim_path, @sprintf("%08d", step))
+    p = Progress(length(1:1:length(groups)) - 1; dt=3.0,
+        desc      = "\e[1;36m[ Info:\e[0m writing",
+        barlen    = 20,
+        barglyphs = BarGlyphs(" ■■  ")
+    )
+    @inbounds Threads.@threads for igroup in eachindex(groups)
+        grp = fid[groups[igroup][1]]
+        coords = permutedims(read(grp["ξ"]))
+        vtk_path = joinpath(anim_path, @sprintf("%08d", igroup))
         vtp_cls  = [MeshCell(PolyData.Verts(), [i]) for i in 1:size(coords, 2)]
         vtk_grid(vtk_path, coords, vtp_cls; compress=true, append=false, ascii=false) do vtk
             for name in keys(grp)
                 name in ("ξ", "time") && continue
-                vtk[name] = read(grp[name])
+                data = read(grp[name])
+                data = ndims(data) == 2 ? permutedims(data) : data
+                vtk[name] = data
             end
-            pvd[t] = vtk # current simulation time
-        end
-
-        t2 = time(); if t2 - t1 > 3.0
-            percen = "animation (.pvd): " * @sprintf("%6.2f%%", (p_iters / p_total) * 100)
-            eta = "eta: " * format_seconds((p_total - p_iters) / (p_iters / (t2 - t_start)))
-            invo_str = "   \e[1;32m⇌\e[0m   "
-            info_con = "\e[1;32m[ Info: \e[0m"
-            print(stdout, "\r\e[2K"); print(stdout, info_con * percen * invo_str * eta)
-            t1 = t2
         end
-        p_iters += 1; p_iters + 1 ≥ p_total && println()
+        next!(p)
     end
-    close(pvd) # close Paraview collection
     _grid2vtr(fid, joinpath(conf.prjdst, "grid.vtr"))
     close(fid) # close HDF5 文件
+    times = map(last, groups)  # 取出排序后的真实时间
+    pvd_path = joinpath(conf.prjdst, "$(conf.prjname).pvd")
+    _write_pvd(pvd_path, times; basedir="animation", ext=".vtp", part=0)
+    return nothing
 end
 
 @views function _grid2vtr(fid, out_file::String)
@@ -97,4 +95,25 @@ end
     end
     vtr["h"] = h
     close(vtr)
+end
+
+function _write_pvd(outpath::AbstractString, times::AbstractVector{<:Real};
+                    basedir::AbstractString = "animation",
+                    ext::AbstractString = ".vtp",
+                    part::Integer = 0)
+    open(outpath, "w") do io
+        println(io, """<?xml version="1.0" encoding="utf-8"?>""")
+        println(io, """<VTKFile type="Collection" version="1.0" byte_order="LittleEndian" compressor="vtkZLibDataCompressor">""")
+        println(io, "  <Collection>")
+        # 用 %.17g 写 timestep，既短又精确（避免 0.14749999999999996 这种串太长）
+        for (i, t) in enumerate(times)
+            step = @sprintf("%08d", i)
+            tstr = @sprintf("%.17g", float(t))
+            fileattr = string(basedir, "/", step, ext)  # 强制使用正斜杠，跨平台稳妥
+            println(io, """    <DataSet timestep="$tstr" part="$part" file="$fileattr"/>""")
+        end
+        println(io, "  </Collection>")
+        println(io, "</VTKFile>")
+    end
+    return outpath
 end

src/precompile.jl

@@ -1,8 +1,8 @@
 @setup_workload begin
     @compile_workload begin
         quiet() do
-            scatter(rand(10, 3))
-            pts2vtp(rand(2, 10); vtp_file=joinpath(tempdir(), "output.vtp"), data=(x=rand(10), y=rand(10)))
+            #scatter(rand(10, 3))
+            # pts2vtp(rand(2, 10); vtp_file=joinpath(tempdir(), "output.vtp"), data=(x=rand(10), y=rand(10)))
         end
     end
 end

src/pts2surf.jl

@@ -62,8 +62,8 @@ function pts2surf(coords::AbstractArray, output_file::String;
     subdomain_num_cubes_per_dim::Int=64,
     output_mesh_smoothing_weights::Bool=true
 )
-    @assert size(coords, 1) ∈ (2, 3) "coords should be a 2D or 3D array"
-    particles = np.array(coords', dtype=np.float64)
+    @assert size(coords, 2) ∈ (2, 3) "coords should be a 2D or 3D array"
+    particles = np.array(coords, dtype=np.float64)
     @assert isfile(output_file) "output_file should be a valid file path"
     file_name = endswith(lowercase(output_file), ".obj") ? output_file : output_file * ".obj"
     mesh_with_data, reconstruction = splashsurf.reconstruction_pipeline(particles,
@@ -177,7 +177,7 @@ function hdf2surf(conf::NamedTuple;
     t1, t_start = time(), time()
     @inbounds for (step, (gname, t)) in enumerate(groups)
         grp = fid[gname]
-        coords = read(grp["ξ"])'  # 读取粒子坐标
+        coords = read(grp["ξ"])  # 读取粒子坐标
         obj_filename = joinpath(obj_path, "$(@sprintf("%08d", step)).obj")
         pts2surf(coords, obj_filename;
             radius=radius,

src/pts2vtp.jl

@@ -19,12 +19,14 @@ Description:
 Generates a `.vtp` file by passing custom fields.
 """
 function pts2vtp(coords; vtp_file="output.vtp", data::NamedTuple=NamedTuple())
-    pts_num = size(coords, 2)
-    size(coords, 1) in [2, 3] || throw(ArgumentError("The input coordinates should be 2D or 3D"))
+    pts_num = size(coords, 1)
+    size(coords, 2) in [2, 3] || throw(ArgumentError("The input coordinates should be 2D or 3D"))
     vtp_cls = [MeshCell(PolyData.Verts(), [i]) for i in 1:pts_num]
-    vtk_grid(vtp_file, coords, vtp_cls; compress=true, append=false, ascii=false) do vtk
+    vtk_grid(vtp_file, coords', vtp_cls; compress=true, append=false, ascii=false) do vtk
         keys(data) ≠ () && for vtp_key in keys(data)
-            vtk[string(vtp_key)] = getfield(data, vtp_key)
+            vdata = getfield(data, vtp_key)
+            vdata = ndims(vdata) == 2 ? permutedims(vdata) : vdata
+            vtk[string(vtp_key)] = vdata
         end
     end
 end