From 1baadc8e1ee32ebd4e5e6a80cada92574f40d35b Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Sun, 18 Jan 2026 23:32:11 +0100 Subject: [PATCH 1/9] Add StillWedge point probe --- README.md | 7 + example/StillWedgeMDBC.jl | 9 +- src/PointMeasurements.jl | 132 +++++++++++++++++ src/ProduceHDFVTK.jl | 193 ++++++++++++++++++------- src/SPHExample.jl | 4 + src/SimulationMetaDataConfiguration.jl | 3 + 6 files changed, 298 insertions(+), 50 deletions(-) create mode 100644 src/PointMeasurements.jl diff --git a/README.md b/README.md index e1bf6bb4..c04a487b 100644 --- a/README.md +++ b/README.md @@ -43,6 +43,7 @@ Example scripts live in `example/` and read geometry from `input/`. The solver c src/ ├── AuxiliaryFunctions.jl # Small helper utilities ├── OpenExternalPrograms.jl # Convenience wrappers for logs and ParaView +├── PointMeasurements.jl # Point probe configuration helpers ├── PreProcess.jl # Load inputs and allocate arrays ├── ProduceHDFVTK.jl # Write simulation data in HDF5/VTK format ├── SPHCellList.jl # Custom neighbour search and time stepping @@ -89,6 +90,12 @@ To color exported cell grids by particle counts, set `ParticleNeighborsPerCell` array that includes each cell's particle count minus one plus the particles in its neighbor stencil. +To record a property at a specific point, populate `PointMeasures` in the +`SimulationMetaData` configuration. Each `PointMeasure` names a location and a +list of variables to interpolate (or fallback to the nearest particle when the +kernel support is empty). These values are written into the main VTKHDF output +as `FieldData` arrays so they travel with each time step. + ## Help Questions or issues can be posted on the GitHub issue tracker. Response times may vary but all feedback is welcome. diff --git a/example/StillWedgeMDBC.jl b/example/StillWedgeMDBC.jl index 76795e98..1833ee1f 100644 --- a/example/StillWedgeMDBC.jl +++ b/example/StillWedgeMDBC.jl @@ -1,4 +1,5 @@ using SPHExample +using StaticArrays let Dimensions = 2 @@ -15,6 +16,13 @@ let VisualizeInParaview=true, ExportSingleVTKHDF=true, ExportGridCells=true, + PointMeasures=[ + PointMeasure( + "Middle", + SVector{Dimensions, FloatType}(0.5, 0.25), + ["Pressure", "Velocity"], + ), + ], OpenLogFile=true, ) @@ -98,4 +106,3 @@ let # display(plt) end - diff --git a/src/PointMeasurements.jl b/src/PointMeasurements.jl new file mode 100644 index 00000000..45444def --- /dev/null +++ b/src/PointMeasurements.jl @@ -0,0 +1,132 @@ +module PointMeasurements + +using LinearAlgebra +using StaticArrays + +using ..SPHKernels: Wᵢⱼ + +export PointMeasure, PointMeasureFieldNames, FillPointMeasureData! + +""" + PointMeasure(Name, Position, Variables) + +Define a probe point at `Position` for measuring `Variables` during output. +""" +struct PointMeasure{Dimensions, FloatType} + Name::String + Position::SVector{Dimensions, FloatType} + Variables::Vector{String} +end + +@inline function To3DPoint(::Val{2}, Point::SVector{2, T}) where {T} + return SVector{3, T}(Point[1], Point[2], zero(T)) +end + +@inline function To3DPoint(::Val{3}, Point::SVector{3, T}) where {T} + return Point +end + +function PointMeasureFieldNames(PointMeasures::Vector{<:PointMeasure}) + FieldNames = String[] + for Measure in PointMeasures + push!(FieldNames, "PointMeasure_$(Measure.Name)_Position") + for Variable in Measure.Variables + push!(FieldNames, "PointMeasure_$(Measure.Name)_$(Variable)") + end + end + return FieldNames +end + +function FillPointMeasureData!(FieldData, PointMeasures::Vector{<:PointMeasure}, + SimParticles, SimKernel, Dimensions; + field_map, vector_fields) + resize!(FieldData, 0) + if isempty(PointMeasures) + return FieldData + end + + for Measure in PointMeasures + AppendPointMeasureResults!(FieldData, Measure, SimParticles, SimKernel, Dimensions; + field_map = field_map, vector_fields = vector_fields) + end + + return FieldData +end + +function AppendPointMeasureResults!(FieldData, Measure::PointMeasure{D, T}, + SimParticles, SimKernel, Dimensions; + field_map, vector_fields) where {D, T} + Positions = SimParticles.Position + NumberOfPoints = length(Positions) + if NumberOfPoints == 0 + error("PointMeasure cannot evaluate without particles.") + end + + Position = Measure.Position + push!(FieldData, To3DPoint(Val(Dimensions), Position)) + + Variables = Measure.Variables + Accumulators = Vector{Any}(undef, length(Variables)) + InterpolateMask = BitVector(undef, length(Variables)) + Sources = Vector{Any}(undef, length(Variables)) + + for (Index, Variable) in pairs(Variables) + Property = get(field_map, Variable, Symbol(Variable)) + if !hasproperty(SimParticles, Property) + error("PointMeasure includes $(Variable) but SimParticles has no field $(Property).") + end + Source = getproperty(SimParticles, Property) + Sources[Index] = Source + if Variable in vector_fields + Accumulators[Index] = zero(eltype(Source)) + InterpolateMask[Index] = eltype(eltype(Source)) <: AbstractFloat + else + Accumulators[Index] = zero(eltype(Source)) + InterpolateMask[Index] = eltype(Source) <: AbstractFloat + end + end + + NearestIndex = 1 + NearestDistance = typemax(T) + WeightSum = zero(T) + + @inbounds for ParticleIndex in eachindex(Positions) + Offset = Position - Positions[ParticleIndex] + Distance = norm(Offset) + if Distance < NearestDistance + NearestDistance = Distance + NearestIndex = ParticleIndex + end + q = clamp(Distance * SimKernel.h⁻¹, zero(T), T(2)) + if q <= T(2) + Weight = Wᵢⱼ(SimKernel, q) + if Weight != zero(T) + WeightSum += Weight + for VariableIndex in eachindex(Variables) + if InterpolateMask[VariableIndex] + Accumulators[VariableIndex] += Sources[VariableIndex][ParticleIndex] * Weight + end + end + end + end + end + + for (Index, Variable) in pairs(Variables) + Source = Sources[Index] + Value = if InterpolateMask[Index] && WeightSum > zero(T) + Accumulators[Index] / WeightSum + else + Source[NearestIndex] + end + + if Variable in vector_fields + push!(FieldData, To3DPoint(Val(Dimensions), Value)) + else + push!(FieldData, Value) + end + end + + return FieldData +end + +end diff --git a/src/ProduceHDFVTK.jl b/src/ProduceHDFVTK.jl index b1be0eb5..575b4e90 100644 --- a/src/ProduceHDFVTK.jl +++ b/src/ProduceHDFVTK.jl @@ -23,20 +23,22 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, using StaticArrays using ..AuxiliaryFunctions: to_3d! + using ..PointMeasurements: PointMeasureFieldNames, FillPointMeasureData! const idType = Int64 const fType = Float64 - struct ParticleSnapshot{P, V} + struct ParticleSnapshot{P, V, F} positions::P output_data::V + field_data::F end - struct ParticleWriteJob{P, V} + struct ParticleWriteJob{P, V, F} iteration::Int time::AbstractFloat - snapshot::ParticleSnapshot{P, V} + snapshot::ParticleSnapshot{P, V, F} end struct GridWriteJob{N} @@ -148,8 +150,10 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, return payload end - function SaveVTKHDF(fid_vector, index, filepath, points, variable_names = String[], args...) + function SaveVTKHDF(fid_vector, index, filepath, points, variable_names = String[], args...; + field_data_names = String[], field_data = Any[]) @assert length(variable_names) == length(args) "Same number of variable_names as args is necessary" + @assert length(field_data_names) == length(field_data) "Same number of field_data_names as field_data is necessary" io = h5open(filepath, "w") gtop = HDF5.create_group(io, "VTKHDF") @@ -168,6 +172,18 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, end end + # Field data + let g = HDF5.create_group(gtop, "FieldData") + for i ∈ eachindex(field_data_names) + value = field_data[i] + if value isa AbstractArray + g[field_data_names[i]] = reinterpret(reshape, eltype(value), value) + else + g[field_data_names[i]] = [value] + end + end + end + # Vertices: 1 point per cell let g = HDF5.create_group(gtop, "Vertices") g["NumberOfCells"] = [np] @@ -196,8 +212,11 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, vtk_file_type = "PolyData", idType = Int64, fType = Float64, - cell_data_names = ["CellData"]) + cell_data_names = ["CellData"], + field_data_names = String[], + field_data_init = Any[]) @assert length(variable_names) == length(args) "Same number of variable_names as args is necessary" + @assert length(field_data_names) == length(field_data_init) "Same number of field_data_names as field_data_init is necessary" # Write version of VTKHDF format as an attribute HDF5.attrs(root)["Version"] = Int32.([2, 3]) @@ -231,6 +250,20 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, HDF5.create_dataset(pData, var_name, arg_val_type, ((0,),(-1,)), chunk=(chunk_size,)) end end + + FieldData = HDF5.create_group(root, "FieldData") + for i ∈ eachindex(field_data_names) + field_name = field_data_names[i] + field_value = field_data_init[i] + field_value_type = field_value isa Number ? typeof(field_value) : eltype(field_value) + field_value_length = field_value isa AbstractArray ? length(field_value) : 1 + + if field_value_length == 3 + HDF5.create_dataset(FieldData, field_name, field_value_type, ((3,0),(3,-1)), chunk=(3,chunk_size)) + else + HDF5.create_dataset(FieldData, field_name, field_value_type, ((0,),(-1,)), chunk=(chunk_size,)) + end + end elseif vtk_file_type == "UnstructuredGrid" HDF5.create_dataset(root, "Connectivity" , idType , ((0,),(-1,)), chunk=(chunk_size,)) HDF5.create_dataset(root, "NumberOfCells" , idType , ((0,),(-1,)), chunk=(chunk_size,)) @@ -238,7 +271,19 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, HDF5.create_dataset(root, "Offsets" , idType , ((0,),(-1,)), chunk=(chunk_size,)) HDF5.create_dataset(root, "Types" , UInt8 , ((0,),(-1,)), chunk=(chunk_size,)) #Must be UInt8 - FieldData = HDF5.create_group(root, "FieldData") #Currently just empty group + FieldData = HDF5.create_group(root, "FieldData") + for i ∈ eachindex(field_data_names) + field_name = field_data_names[i] + field_value = field_data_init[i] + field_value_type = field_value isa Number ? typeof(field_value) : eltype(field_value) + field_value_length = field_value isa AbstractArray ? length(field_value) : 1 + + if field_value_length == 3 + HDF5.create_dataset(FieldData, field_name, field_value_type, ((3,0),(3,-1)), chunk=(3,chunk_size)) + else + HDF5.create_dataset(FieldData, field_name, field_value_type, ((0,),(-1,)), chunk=(chunk_size,)) + end + end CellData = HDF5.create_group(root, "CellData") for name in cell_data_names @@ -254,7 +299,8 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, function GenerateStepStructure(root, variable_names = String[], args...; vtk_file_type = "PolyData", chunk_size = 1000, - cell_data_names = ["CellData"]) + cell_data_names = ["CellData"], + field_data_names = String[]) steps = HDF5.create_group(root, "Steps") NSteps, _ = HDF5.create_attribute(steps, "NSteps", Int32) @@ -271,6 +317,10 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, for name in nTopoDSs HDF5.create_dataset(steps, name, idType, ((4,0),(4, -1)), chunk=(4, chunk_size)) end + fData = HDF5.create_group(steps, "FieldDataOffsets") + for name in field_data_names + HDF5.create_dataset(fData, name, idType, ((0,),(-1,)), chunk=(chunk_size,)) + end elseif vtk_file_type == "UnstructuredGrid" nTopoDSs = ["CellOffsets", "ConnectivityIdOffsets"] for name in nTopoDSs @@ -281,6 +331,10 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, HDF5.create_dataset(cData, name, idType, ((0,),(-1,)), chunk=(chunk_size,)) end + fData = HDF5.create_group(steps, "FieldDataOffsets") + for name in field_data_names + HDF5.create_dataset(fData, name, idType, ((0,),(-1,)), chunk=(chunk_size,)) + end end pData = HDF5.create_group(steps, "PointDataOffsets") @@ -294,7 +348,8 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, end - function AppendVTKHDFData(root, newStep, Positions, variable_names, args...) + function AppendVTKHDFData(root, newStep, Positions, variable_names, args...; + field_data_names = String[], field_data = Any[]) steps = root["Steps"] # To update attributes, this is the best way I've found so far @@ -359,6 +414,24 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, root["PointData"][var_name][PointsStartIndex:(PointsStartIndex + PositionLength - 1)] = arg end end + + if !isempty(field_data_names) + for i ∈ eachindex(field_data_names) + field_name = field_data_names[i] + field_value = field_data[i] + if field_value isa AbstractArray + FieldDataStartIndex = size(root["FieldData"][field_name], 2) + 1 + HDF5.set_extent_dims(root["FieldData"][field_name], (length(field_value), FieldDataStartIndex)) + root["FieldData"][field_name][:, FieldDataStartIndex] = field_value + else + FieldDataStartIndex = length(root["FieldData"][field_name]) + 1 + HDF5.set_extent_dims(root["FieldData"][field_name], (FieldDataStartIndex,)) + root["FieldData"][field_name][FieldDataStartIndex] = field_value + end + HDF5.set_extent_dims(steps["FieldDataOffsets"][field_name], (length(steps["FieldDataOffsets"][field_name]) + 1,)) + steps["FieldDataOffsets"][field_name][end] = Int(FieldDataStartIndex - 1) + end + end connectivities = ["Vertices", "Lines", "Polygons", "Strips"] @@ -548,6 +621,29 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, grid_filename = (iter) -> "$(grid_savepath)_$(lpad(iter,6,"0")).vtkhdf" output_vars = SimMetaData.OutputVariables + point_measure_names = PointMeasureFieldNames(SimMetaData.PointMeasures) + + vector_fields = Set([ + "KernelGradient", + "Velocity", + "Acceleration", + "GhostPoints", + "GhostNormals", + ]) + field_map = Dict( + "Kernel" => :Kernel, + "KernelGradient" => :KernelGradient, + "Density" => :Density, + "Pressure" => :Pressure, + "Velocity" => :Velocity, + "Acceleration" => :Acceleration, + "BoundaryBool" => :BoundaryBool, + "ID" => :ID, + "Type" => :Type, + "GroupMarker" => :GroupMarker, + "GhostPoints" => :GhostPoints, + "GhostNormals" => :GhostNormals, + ) # Initialize storage for file handles file_handles = if !SimMetaData.ExportSingleVTKHDF @@ -566,20 +662,6 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, OutputVTKHDF = h5open("$(particle_savepath).vtkhdf", "w") root = HDF5.create_group(OutputVTKHDF, "VTKHDF") - field_map = Dict( - "Kernel" => :Kernel, - "KernelGradient" => :KernelGradient, - "Density" => :Density, - "Pressure" => :Pressure, - "Velocity" => :Velocity, - "Acceleration" => :Acceleration, - "BoundaryBool" => :BoundaryBool, - "ID" => :ID, - "Type" => :Type, - "GroupMarker" => :GroupMarker, - "GhostPoints" => :GhostPoints, - "GhostNormals" => :GhostNormals, - ) output_data_init = Vector{Any}(undef, length(output_vars)) for (i, name) in pairs(output_vars) prop = get(field_map, name, nothing) @@ -593,8 +675,29 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, end end - GenerateGeometryStructure(root, output_vars, output_data_init...; chunk_size=1024) - GenerateStepStructure(root, output_vars, output_data_init...) + point_measure_init = FillPointMeasureData!( + Any[], + SimMetaData.PointMeasures, + SimParticles, + SimKernel, + Dimensions; + field_map = field_map, + vector_fields = vector_fields, + ) + GenerateGeometryStructure( + root, + output_vars, + output_data_init...; + chunk_size = 1024, + field_data_names = point_measure_names, + field_data_init = point_measure_init, + ) + GenerateStepStructure( + root, + output_vars, + output_data_init...; + field_data_names = point_measure_names, + ) # Initialize grid file if needed if SimMetaData.ExportGridCells @@ -622,28 +725,6 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, end end - vector_fields = Set([ - "KernelGradient", - "Velocity", - "Acceleration", - "GhostPoints", - "GhostNormals", - ]) - field_map = Dict( - "Kernel" => :Kernel, - "KernelGradient" => :KernelGradient, - "Density" => :Density, - "Pressure" => :Pressure, - "Velocity" => :Velocity, - "Acceleration" => :Acceleration, - "BoundaryBool" => :BoundaryBool, - "ID" => :ID, - "Type" => :Type, - "GroupMarker" => :GroupMarker, - "GhostPoints" => :GhostPoints, - "GhostNormals" => :GhostNormals, - ) - T = eltype(eltype(SimParticles.Position)) function allocate_particle_snapshot() n = length(SimParticles.Position) @@ -663,7 +744,8 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, output_data[i] = similar(src, n) end end - return ParticleSnapshot(positions, output_data) + field_data = Any[] + return ParticleSnapshot(positions, output_data, field_data) end function fill_particle_snapshot!(snapshot) @@ -700,6 +782,15 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, copy!(buf, src) end end + FillPointMeasureData!( + snapshot.field_data, + SimMetaData.PointMeasures, + SimParticles, + SimKernel, + Dimensions; + field_map = field_map, + vector_fields = vector_fields, + ) return snapshot end @@ -719,7 +810,9 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, particle_filename(job.iteration), snapshot.positions, output_vars, - snapshot.output_data..., + snapshot.output_data...; + field_data_names = point_measure_names, + field_data = snapshot.field_data, ) else AppendVTKHDFData( @@ -727,7 +820,9 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, job.time, snapshot.positions, output_vars, - snapshot.output_data..., + snapshot.output_data...; + field_data_names = point_measure_names, + field_data = snapshot.field_data, ) end put!(buffer_pool, snapshot) diff --git a/src/SPHExample.jl b/src/SPHExample.jl index 28ef3ae7..a8817a4d 100644 --- a/src/SPHExample.jl +++ b/src/SPHExample.jl @@ -3,6 +3,7 @@ module SPHExample include("AuxiliaryFunctions.jl"); include("SPHKernels.jl") include("SPHViscosityModels.jl") + include("PointMeasurements.jl") include("ProduceHDFVTK.jl") include("SimulationGeometry.jl") include("SimulationMetaDataConfiguration.jl"); @@ -23,6 +24,9 @@ module SPHExample using .SPHKernels export SPHKernel, SPHKernelInstance, WendlandC2, CubicSpline, Wᵢⱼ, ∇Wᵢⱼ, tensile_correction + using .PointMeasurements + export PointMeasure, PointMeasureFieldNames, FillPointMeasureData! + using .SPHViscosityModels export SPHViscosity, ZeroViscosity, ArtificialViscosity, Laminar, LaminarSPS, compute_viscosity diff --git a/src/SimulationMetaDataConfiguration.jl b/src/SimulationMetaDataConfiguration.jl index c0ce65df..437950b2 100644 --- a/src/SimulationMetaDataConfiguration.jl +++ b/src/SimulationMetaDataConfiguration.jl @@ -3,6 +3,8 @@ module SimulationMetaDataConfiguration using Parameters using TimerOutputs +using ..PointMeasurements: PointMeasure + export SimulationMetaData, UpdateMetaData!, ShiftingMode, NoShifting, PlanarShifting, KernelOutputMode, NoKernelOutput, StoreKernelOutput, @@ -56,6 +58,7 @@ struct StoreLog <: LogMode end "Type", "GroupMarker", ] + PointMeasures::Vector{PointMeasure} = PointMeasure[] OpenLogFile::Bool = true Δx::FloatType = zero(FloatType) end From 6fd0a4684a99a920306144edc91e2579d9060266 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Mon, 19 Jan 2026 00:02:50 +0100 Subject: [PATCH 2/9] Use neighbor data for point measures --- src/PointMeasurements.jl | 90 ++++++++++++++++++++++++++++++++-------- src/ProduceHDFVTK.jl | 7 ++-- src/SPHCellList.jl | 9 +++- src/SPHExample.jl | 2 +- 4 files changed, 86 insertions(+), 22 deletions(-) diff --git a/src/PointMeasurements.jl b/src/PointMeasurements.jl index 45444def..addb1f46 100644 --- a/src/PointMeasurements.jl +++ b/src/PointMeasurements.jl @@ -4,6 +4,7 @@ using LinearAlgebra using StaticArrays using ..SPHKernels: Wᵢⱼ +using ..SPHNeighborList: MapFloor export PointMeasure, PointMeasureFieldNames, FillPointMeasureData! @@ -39,7 +40,8 @@ end function FillPointMeasureData!(FieldData, PointMeasures::Vector{<:PointMeasure}, SimParticles, SimKernel, Dimensions; - field_map, vector_fields) + field_map, vector_fields, + neighbor_data = nothing) resize!(FieldData, 0) if isempty(PointMeasures) return FieldData @@ -47,15 +49,31 @@ function FillPointMeasureData!(FieldData, PointMeasures::Vector{<:PointMeasure}, for Measure in PointMeasures AppendPointMeasureResults!(FieldData, Measure, SimParticles, SimKernel, Dimensions; - field_map = field_map, vector_fields = vector_fields) + field_map = field_map, + vector_fields = vector_fields, + neighbor_data = neighbor_data) end return FieldData end +function FindNearestIndex(Positions, Position) + NearestIndex = 1 + NearestDistance = typemax(eltype(Position)) + @inbounds for ParticleIndex in eachindex(Positions) + Distance = norm(Position - Positions[ParticleIndex]) + if Distance < NearestDistance + NearestDistance = Distance + NearestIndex = ParticleIndex + end + end + return NearestIndex, NearestDistance +end + function AppendPointMeasureResults!(FieldData, Measure::PointMeasure{D, T}, SimParticles, SimKernel, Dimensions; - field_map, vector_fields) where {D, T} + field_map, vector_fields, + neighbor_data) where {D, T} Positions = SimParticles.Position NumberOfPoints = length(Positions) if NumberOfPoints == 0 @@ -90,21 +108,59 @@ function AppendPointMeasureResults!(FieldData, Measure::PointMeasure{D, T}, NearestDistance = typemax(T) WeightSum = zero(T) - @inbounds for ParticleIndex in eachindex(Positions) - Offset = Position - Positions[ParticleIndex] - Distance = norm(Offset) - if Distance < NearestDistance - NearestDistance = Distance - NearestIndex = ParticleIndex + HasCandidates = false + if neighbor_data !== nothing + cell_dict = neighbor_data.cell_dict + particle_ranges = neighbor_data.particle_ranges + full_stencil = neighbor_data.full_stencil + cell_index = CartesianIndex(map(x -> MapFloor(x, SimKernel.H⁻¹), Tuple(Position))) + + @inbounds for offset in full_stencil + neighbor_cell = cell_index + offset + neighbor_index = get(cell_dict, neighbor_cell, 0) + if neighbor_index != 0 + start_index = particle_ranges[neighbor_index] + end_index = particle_ranges[neighbor_index + 1] - 1 + if start_index <= end_index + HasCandidates = true + for particle_index in start_index:end_index + Offset = Position - Positions[particle_index] + Distance = norm(Offset) + if Distance < NearestDistance + NearestDistance = Distance + NearestIndex = particle_index + end + q = clamp(Distance * SimKernel.h⁻¹, zero(T), T(2)) + if q <= T(2) + Weight = Wᵢⱼ(SimKernel, q) + if Weight != zero(T) + WeightSum += Weight + for VariableIndex in eachindex(Variables) + if InterpolateMask[VariableIndex] + Accumulators[VariableIndex] += Sources[VariableIndex][particle_index] * Weight + end + end + end + end + end + end + end end - q = clamp(Distance * SimKernel.h⁻¹, zero(T), T(2)) - if q <= T(2) - Weight = Wᵢⱼ(SimKernel, q) - if Weight != zero(T) - WeightSum += Weight - for VariableIndex in eachindex(Variables) - if InterpolateMask[VariableIndex] - Accumulators[VariableIndex] += Sources[VariableIndex][ParticleIndex] * Weight + end + + if !HasCandidates + NearestIndex, NearestDistance = FindNearestIndex(Positions, Position) + @inbounds for ParticleIndex in eachindex(Positions) + Distance = norm(Position - Positions[ParticleIndex]) + q = clamp(Distance * SimKernel.h⁻¹, zero(T), T(2)) + if q <= T(2) + Weight = Wᵢⱼ(SimKernel, q) + if Weight != zero(T) + WeightSum += Weight + for VariableIndex in eachindex(Variables) + if InterpolateMask[VariableIndex] + Accumulators[VariableIndex] += Sources[VariableIndex][ParticleIndex] * Weight + end end end end diff --git a/src/ProduceHDFVTK.jl b/src/ProduceHDFVTK.jl index 575b4e90..2da5022e 100644 --- a/src/ProduceHDFVTK.jl +++ b/src/ProduceHDFVTK.jl @@ -748,7 +748,7 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, return ParticleSnapshot(positions, output_data, field_data) end - function fill_particle_snapshot!(snapshot) + function fill_particle_snapshot!(snapshot; neighbor_data = nothing) if Dimensions == 2 to_3d!(snapshot.positions, SimParticles.Position) else @@ -790,6 +790,7 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, Dimensions; field_map = field_map, vector_fields = vector_fields, + neighbor_data = neighbor_data, ) return snapshot end @@ -849,9 +850,9 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, end end - function enqueue_particle_data(iteration) + function enqueue_particle_data(iteration; neighbor_data = nothing) snapshot = take!(buffer_pool) - fill_particle_snapshot!(snapshot) + fill_particle_snapshot!(snapshot; neighbor_data = neighbor_data) job = ParticleWriteJob(iteration, SimMetaData.TotalTime, snapshot) put!(job_channel, job) end diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index eb55a820..b6504538 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -884,7 +884,14 @@ using LinearAlgebra end @timeit SimMetaData.HourGlass "13 Save Particle Data" begin - output.enqueue_particles(SimMetaData.OutputIterationCounter) + output.enqueue_particles( + SimMetaData.OutputIterationCounter; + neighbor_data = ( + particle_ranges = ParticleRanges, + cell_dict = CellDict, + full_stencil = FullStencil, + ), + ) output.enqueue_grid(SimMetaData.OutputIterationCounter, UniqueCellsView, cell_particle_counts=cell_particle_counts, cell_neighbor_counts=cell_neighbor_counts) end diff --git a/src/SPHExample.jl b/src/SPHExample.jl index a8817a4d..2bf4a2c2 100644 --- a/src/SPHExample.jl +++ b/src/SPHExample.jl @@ -3,6 +3,7 @@ module SPHExample include("AuxiliaryFunctions.jl"); include("SPHKernels.jl") include("SPHViscosityModels.jl") + include("SPHNeighborList.jl") include("PointMeasurements.jl") include("ProduceHDFVTK.jl") include("SimulationGeometry.jl") @@ -14,7 +15,6 @@ module SPHExample include("PreProcess.jl"); include("OpenExternalPrograms.jl") include("SPHDensityDiffusionModels.jl") - include("SPHNeighborList.jl") include("SPHCellList.jl") #Must be last # Re-export desired functions from each submodule From feb42c68992dea5116cee34d9ba2a0b504d997ec Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Mon, 19 Jan 2026 00:02:55 +0100 Subject: [PATCH 3/9] Remove brute-force probe fallback --- src/PointMeasurements.jl | 91 ++++++++++++++-------------------------- src/SPHCellList.jl | 9 +++- 2 files changed, 40 insertions(+), 60 deletions(-) diff --git a/src/PointMeasurements.jl b/src/PointMeasurements.jl index addb1f46..43feaff4 100644 --- a/src/PointMeasurements.jl +++ b/src/PointMeasurements.jl @@ -46,6 +46,9 @@ function FillPointMeasureData!(FieldData, PointMeasures::Vector{<:PointMeasure}, if isempty(PointMeasures) return FieldData end + if neighbor_data === nothing + error("PointMeasure evaluation requires neighbor data; pass `neighbor_data` from the cell list.") + end for Measure in PointMeasures AppendPointMeasureResults!(FieldData, Measure, SimParticles, SimKernel, Dimensions; @@ -57,19 +60,6 @@ function FillPointMeasureData!(FieldData, PointMeasures::Vector{<:PointMeasure}, return FieldData end -function FindNearestIndex(Positions, Position) - NearestIndex = 1 - NearestDistance = typemax(eltype(Position)) - @inbounds for ParticleIndex in eachindex(Positions) - Distance = norm(Position - Positions[ParticleIndex]) - if Distance < NearestDistance - NearestDistance = Distance - NearestIndex = ParticleIndex - end - end - return NearestIndex, NearestDistance -end - function AppendPointMeasureResults!(FieldData, Measure::PointMeasure{D, T}, SimParticles, SimKernel, Dimensions; field_map, vector_fields, @@ -108,37 +98,35 @@ function AppendPointMeasureResults!(FieldData, Measure::PointMeasure{D, T}, NearestDistance = typemax(T) WeightSum = zero(T) + cell_dict = neighbor_data.cell_dict + particle_ranges = neighbor_data.particle_ranges + full_stencil = neighbor_data.full_stencil + cell_index = CartesianIndex(map(x -> MapFloor(x, SimKernel.H⁻¹), Tuple(Position))) + HasCandidates = false - if neighbor_data !== nothing - cell_dict = neighbor_data.cell_dict - particle_ranges = neighbor_data.particle_ranges - full_stencil = neighbor_data.full_stencil - cell_index = CartesianIndex(map(x -> MapFloor(x, SimKernel.H⁻¹), Tuple(Position))) - - @inbounds for offset in full_stencil - neighbor_cell = cell_index + offset - neighbor_index = get(cell_dict, neighbor_cell, 0) - if neighbor_index != 0 - start_index = particle_ranges[neighbor_index] - end_index = particle_ranges[neighbor_index + 1] - 1 - if start_index <= end_index - HasCandidates = true - for particle_index in start_index:end_index - Offset = Position - Positions[particle_index] - Distance = norm(Offset) - if Distance < NearestDistance - NearestDistance = Distance - NearestIndex = particle_index - end - q = clamp(Distance * SimKernel.h⁻¹, zero(T), T(2)) - if q <= T(2) - Weight = Wᵢⱼ(SimKernel, q) - if Weight != zero(T) - WeightSum += Weight - for VariableIndex in eachindex(Variables) - if InterpolateMask[VariableIndex] - Accumulators[VariableIndex] += Sources[VariableIndex][particle_index] * Weight - end + @inbounds for offset in full_stencil + neighbor_cell = cell_index + offset + neighbor_index = get(cell_dict, neighbor_cell, 0) + if neighbor_index != 0 + start_index = particle_ranges[neighbor_index] + end_index = particle_ranges[neighbor_index + 1] - 1 + if start_index <= end_index + HasCandidates = true + for particle_index in start_index:end_index + Offset = Position - Positions[particle_index] + Distance = norm(Offset) + if Distance < NearestDistance + NearestDistance = Distance + NearestIndex = particle_index + end + q = clamp(Distance * SimKernel.h⁻¹, zero(T), T(2)) + if q <= T(2) + Weight = Wᵢⱼ(SimKernel, q) + if Weight != zero(T) + WeightSum += Weight + for VariableIndex in eachindex(Variables) + if InterpolateMask[VariableIndex] + Accumulators[VariableIndex] += Sources[VariableIndex][particle_index] * Weight end end end @@ -149,22 +137,7 @@ function AppendPointMeasureResults!(FieldData, Measure::PointMeasure{D, T}, end if !HasCandidates - NearestIndex, NearestDistance = FindNearestIndex(Positions, Position) - @inbounds for ParticleIndex in eachindex(Positions) - Distance = norm(Position - Positions[ParticleIndex]) - q = clamp(Distance * SimKernel.h⁻¹, zero(T), T(2)) - if q <= T(2) - Weight = Wᵢⱼ(SimKernel, q) - if Weight != zero(T) - WeightSum += Weight - for VariableIndex in eachindex(Variables) - if InterpolateMask[VariableIndex] - Accumulators[VariableIndex] += Sources[VariableIndex][ParticleIndex] * Weight - end - end - end - end - end + error("PointMeasure found no neighbor candidates; adjust probe position or cell list.") end for (Index, Variable) in pairs(Variables) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index b6504538..094a3076 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -826,7 +826,14 @@ using LinearAlgebra # Save initial state, use 1 else this cannot be used to index fid vector SimMetaData.OutputIterationCounter = 1 - output.enqueue_particles(SimMetaData.OutputIterationCounter) + output.enqueue_particles( + SimMetaData.OutputIterationCounter; + neighbor_data = ( + particle_ranges = ParticleRanges, + cell_dict = CellDict, + full_stencil = FullStencil, + ), + ) if SimMetaData.IndexCounter > 0 unique_cells_view = view(UniqueCells, 1:SimMetaData.IndexCounter) cell_particle_counts = nothing From 1b42cdb91d7c7cceb3e41e3244c4c39c4b8d2d06 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Mon, 19 Jan 2026 00:03:03 +0100 Subject: [PATCH 4/9] Output point measures as polydata --- README.md | 6 +- src/PointMeasurements.jl | 58 ++++++++++- src/ProduceHDFVTK.jl | 214 ++++++++++++++++++++------------------- 3 files changed, 171 insertions(+), 107 deletions(-) diff --git a/README.md b/README.md index c04a487b..4d678b3b 100644 --- a/README.md +++ b/README.md @@ -92,9 +92,9 @@ one plus the particles in its neighbor stencil. To record a property at a specific point, populate `PointMeasures` in the `SimulationMetaData` configuration. Each `PointMeasure` names a location and a -list of variables to interpolate (or fallback to the nearest particle when the -kernel support is empty). These values are written into the main VTKHDF output -as `FieldData` arrays so they travel with each time step. +list of variables to interpolate within the neighbor stencil. These values are +written into a dedicated `*_PointMeasures.vtkhdf` file so you can load probe +points alongside the particle data. ## Help diff --git a/src/PointMeasurements.jl b/src/PointMeasurements.jl index 43feaff4..ee06f1f9 100644 --- a/src/PointMeasurements.jl +++ b/src/PointMeasurements.jl @@ -6,7 +6,8 @@ using StaticArrays using ..SPHKernels: Wᵢⱼ using ..SPHNeighborList: MapFloor -export PointMeasure, PointMeasureFieldNames, FillPointMeasureData! +export PointMeasure, PointMeasureFieldNames, PointMeasureOutputVariableNames, + FillPointMeasureData!, FillPointMeasureSnapshot! """ PointMeasure(Name, Position, Variables) @@ -38,6 +39,20 @@ function PointMeasureFieldNames(PointMeasures::Vector{<:PointMeasure}) return FieldNames end +function PointMeasureOutputVariableNames(PointMeasures::Vector{<:PointMeasure}) + OutputNames = String[] + Seen = Set{String}() + for Measure in PointMeasures + for Variable in Measure.Variables + if !(Variable in Seen) + push!(OutputNames, Variable) + push!(Seen, Variable) + end + end + end + return OutputNames +end + function FillPointMeasureData!(FieldData, PointMeasures::Vector{<:PointMeasure}, SimParticles, SimKernel, Dimensions; field_map, vector_fields, @@ -158,4 +173,45 @@ function AppendPointMeasureResults!(FieldData, Measure::PointMeasure{D, T}, return FieldData end +function FillPointMeasureSnapshot!(snapshot_positions, snapshot_output_data, + PointMeasures::Vector{<:PointMeasure}, + OutputVariables::Vector{String}, + SimParticles, SimKernel, Dimensions; + field_map, vector_fields, + neighbor_data = nothing) + if isempty(PointMeasures) + return snapshot_positions, snapshot_output_data + end + if neighbor_data === nothing + error("PointMeasure evaluation requires neighbor data; pass `neighbor_data` from the cell list.") + end + + resize!(snapshot_positions, length(PointMeasures)) + for Output in snapshot_output_data + resize!(Output, length(PointMeasures)) + fill!(Output, zero(eltype(Output))) + end + + output_index = Dict(name => idx for (idx, name) in pairs(OutputVariables)) + + for (MeasureIndex, Measure) in pairs(PointMeasures) + snapshot_positions[MeasureIndex] = To3DPoint(Val(Dimensions), Measure.Position) + + buffer = Any[] + AppendPointMeasureResults!(buffer, Measure, SimParticles, SimKernel, Dimensions; + field_map = field_map, + vector_fields = vector_fields, + neighbor_data = neighbor_data) + + data_index = 2 + for Variable in Measure.Variables + index = output_index[Variable] + snapshot_output_data[index][MeasureIndex] = buffer[data_index] + data_index += 1 + end + end + + return snapshot_positions, snapshot_output_data +end + end diff --git a/src/ProduceHDFVTK.jl b/src/ProduceHDFVTK.jl index 2da5022e..bb17b1f1 100644 --- a/src/ProduceHDFVTK.jl +++ b/src/ProduceHDFVTK.jl @@ -23,22 +23,27 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, using StaticArrays using ..AuxiliaryFunctions: to_3d! - using ..PointMeasurements: PointMeasureFieldNames, FillPointMeasureData! + using ..PointMeasurements: PointMeasureOutputVariableNames, FillPointMeasureSnapshot! const idType = Int64 const fType = Float64 - struct ParticleSnapshot{P, V, F} + struct ParticleSnapshot{P, V, M} positions::P output_data::V - field_data::F + point_measure_snapshot::M end - struct ParticleWriteJob{P, V, F} + struct ParticleWriteJob{P, V, M} iteration::Int time::AbstractFloat - snapshot::ParticleSnapshot{P, V, F} + snapshot::ParticleSnapshot{P, V, M} + end + + struct PointMeasureSnapshot{P, V} + positions::P + output_data::V end struct GridWriteJob{N} @@ -150,10 +155,8 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, return payload end - function SaveVTKHDF(fid_vector, index, filepath, points, variable_names = String[], args...; - field_data_names = String[], field_data = Any[]) + function SaveVTKHDF(fid_vector, index, filepath, points, variable_names = String[], args...) @assert length(variable_names) == length(args) "Same number of variable_names as args is necessary" - @assert length(field_data_names) == length(field_data) "Same number of field_data_names as field_data is necessary" io = h5open(filepath, "w") gtop = HDF5.create_group(io, "VTKHDF") @@ -172,18 +175,6 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, end end - # Field data - let g = HDF5.create_group(gtop, "FieldData") - for i ∈ eachindex(field_data_names) - value = field_data[i] - if value isa AbstractArray - g[field_data_names[i]] = reinterpret(reshape, eltype(value), value) - else - g[field_data_names[i]] = [value] - end - end - end - # Vertices: 1 point per cell let g = HDF5.create_group(gtop, "Vertices") g["NumberOfCells"] = [np] @@ -212,11 +203,8 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, vtk_file_type = "PolyData", idType = Int64, fType = Float64, - cell_data_names = ["CellData"], - field_data_names = String[], - field_data_init = Any[]) + cell_data_names = ["CellData"]) @assert length(variable_names) == length(args) "Same number of variable_names as args is necessary" - @assert length(field_data_names) == length(field_data_init) "Same number of field_data_names as field_data_init is necessary" # Write version of VTKHDF format as an attribute HDF5.attrs(root)["Version"] = Int32.([2, 3]) @@ -252,18 +240,6 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, end FieldData = HDF5.create_group(root, "FieldData") - for i ∈ eachindex(field_data_names) - field_name = field_data_names[i] - field_value = field_data_init[i] - field_value_type = field_value isa Number ? typeof(field_value) : eltype(field_value) - field_value_length = field_value isa AbstractArray ? length(field_value) : 1 - - if field_value_length == 3 - HDF5.create_dataset(FieldData, field_name, field_value_type, ((3,0),(3,-1)), chunk=(3,chunk_size)) - else - HDF5.create_dataset(FieldData, field_name, field_value_type, ((0,),(-1,)), chunk=(chunk_size,)) - end - end elseif vtk_file_type == "UnstructuredGrid" HDF5.create_dataset(root, "Connectivity" , idType , ((0,),(-1,)), chunk=(chunk_size,)) HDF5.create_dataset(root, "NumberOfCells" , idType , ((0,),(-1,)), chunk=(chunk_size,)) @@ -272,18 +248,6 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, HDF5.create_dataset(root, "Types" , UInt8 , ((0,),(-1,)), chunk=(chunk_size,)) #Must be UInt8 FieldData = HDF5.create_group(root, "FieldData") - for i ∈ eachindex(field_data_names) - field_name = field_data_names[i] - field_value = field_data_init[i] - field_value_type = field_value isa Number ? typeof(field_value) : eltype(field_value) - field_value_length = field_value isa AbstractArray ? length(field_value) : 1 - - if field_value_length == 3 - HDF5.create_dataset(FieldData, field_name, field_value_type, ((3,0),(3,-1)), chunk=(3,chunk_size)) - else - HDF5.create_dataset(FieldData, field_name, field_value_type, ((0,),(-1,)), chunk=(chunk_size,)) - end - end CellData = HDF5.create_group(root, "CellData") for name in cell_data_names @@ -299,8 +263,7 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, function GenerateStepStructure(root, variable_names = String[], args...; vtk_file_type = "PolyData", chunk_size = 1000, - cell_data_names = ["CellData"], - field_data_names = String[]) + cell_data_names = ["CellData"]) steps = HDF5.create_group(root, "Steps") NSteps, _ = HDF5.create_attribute(steps, "NSteps", Int32) @@ -317,10 +280,6 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, for name in nTopoDSs HDF5.create_dataset(steps, name, idType, ((4,0),(4, -1)), chunk=(4, chunk_size)) end - fData = HDF5.create_group(steps, "FieldDataOffsets") - for name in field_data_names - HDF5.create_dataset(fData, name, idType, ((0,),(-1,)), chunk=(chunk_size,)) - end elseif vtk_file_type == "UnstructuredGrid" nTopoDSs = ["CellOffsets", "ConnectivityIdOffsets"] for name in nTopoDSs @@ -331,10 +290,6 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, HDF5.create_dataset(cData, name, idType, ((0,),(-1,)), chunk=(chunk_size,)) end - fData = HDF5.create_group(steps, "FieldDataOffsets") - for name in field_data_names - HDF5.create_dataset(fData, name, idType, ((0,),(-1,)), chunk=(chunk_size,)) - end end pData = HDF5.create_group(steps, "PointDataOffsets") @@ -348,8 +303,7 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, end - function AppendVTKHDFData(root, newStep, Positions, variable_names, args...; - field_data_names = String[], field_data = Any[]) + function AppendVTKHDFData(root, newStep, Positions, variable_names, args...) steps = root["Steps"] # To update attributes, this is the best way I've found so far @@ -415,25 +369,6 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, end end - if !isempty(field_data_names) - for i ∈ eachindex(field_data_names) - field_name = field_data_names[i] - field_value = field_data[i] - if field_value isa AbstractArray - FieldDataStartIndex = size(root["FieldData"][field_name], 2) + 1 - HDF5.set_extent_dims(root["FieldData"][field_name], (length(field_value), FieldDataStartIndex)) - root["FieldData"][field_name][:, FieldDataStartIndex] = field_value - else - FieldDataStartIndex = length(root["FieldData"][field_name]) + 1 - HDF5.set_extent_dims(root["FieldData"][field_name], (FieldDataStartIndex,)) - root["FieldData"][field_name][FieldDataStartIndex] = field_value - end - HDF5.set_extent_dims(steps["FieldDataOffsets"][field_name], (length(steps["FieldDataOffsets"][field_name]) + 1,)) - steps["FieldDataOffsets"][field_name][end] = Int(FieldDataStartIndex - 1) - end - end - - connectivities = ["Vertices", "Lines", "Polygons", "Strips"] for connect in connectivities for dataset in ["NumberOfCells", "NumberOfConnectivityIds", "Offsets", "Connectivity"] @@ -619,9 +554,10 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, # File naming functions particle_filename = (iter) -> "$(particle_savepath)_$(lpad(iter,6,"0")).vtkhdf" grid_filename = (iter) -> "$(grid_savepath)_$(lpad(iter,6,"0")).vtkhdf" + point_filename = (iter) -> "$(particle_savepath)_PointMeasures_$(lpad(iter,6,"0")).vtkhdf" output_vars = SimMetaData.OutputVariables - point_measure_names = PointMeasureFieldNames(SimMetaData.PointMeasures) + point_measure_vars = PointMeasureOutputVariableNames(SimMetaData.PointMeasures) vector_fields = Set([ "KernelGradient", @@ -644,6 +580,13 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, "GhostPoints" => :GhostPoints, "GhostNormals" => :GhostNormals, ) + + T = eltype(eltype(SimParticles.Position)) + to_point3d = if Dimensions == 2 + (point, ::Type{T}) -> SVector{3, T}(point[1], point[2], zero(T)) + else + (point, ::Type{T}) -> point + end # Initialize storage for file handles file_handles = if !SimMetaData.ExportSingleVTKHDF @@ -656,6 +599,7 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, ( particle_files = Vector{HDF5.File}(undef, n_outputs), grid_files = nothing, + point_files = isempty(point_measure_vars) ? nothing : Vector{HDF5.File}(undef, n_outputs), ) else # Single-file mode: handles for both files @@ -675,28 +619,16 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, end end - point_measure_init = FillPointMeasureData!( - Any[], - SimMetaData.PointMeasures, - SimParticles, - SimKernel, - Dimensions; - field_map = field_map, - vector_fields = vector_fields, - ) GenerateGeometryStructure( root, output_vars, output_data_init...; chunk_size = 1024, - field_data_names = point_measure_names, - field_data_init = point_measure_init, ) GenerateStepStructure( root, output_vars, output_data_init...; - field_data_names = point_measure_names, ) # Initialize grid file if needed @@ -719,13 +651,52 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, cell_data_names=cell_data_names, ) - (particle_files = OutputVTKHDF, grid_files = OutputVTKHDFGrid) + point_handle = nothing + if !isempty(point_measure_vars) + OutputVTKHDFPoints = h5open("$(particle_savepath)_PointMeasures.vtkhdf", "w") + point_root = HDF5.create_group(OutputVTKHDFPoints, "VTKHDF") + point_positions = [to_point3d(measure.Position, T) for measure in SimMetaData.PointMeasures] + point_output_init = Vector{Any}(undef, length(point_measure_vars)) + for (index, name) in pairs(point_measure_vars) + if name in vector_fields + point_output_init[index] = fill(zero(SVector{3, T}), length(SimMetaData.PointMeasures)) + else + prop = field_map[name] + src = getproperty(SimParticles, prop) + point_output_init[index] = fill(zero(eltype(src)), length(SimMetaData.PointMeasures)) + end + end + GenerateGeometryStructure(point_root, point_measure_vars, point_output_init...; chunk_size = 1024) + GenerateStepStructure(point_root, point_measure_vars, point_output_init...) + point_handle = (file = OutputVTKHDFPoints, root = point_root) + end + + (particle_files = OutputVTKHDF, grid_files = OutputVTKHDFGrid, point_files = point_handle) else - (particle_files = OutputVTKHDF, grid_files = nothing) + point_handle = nothing + if !isempty(point_measure_vars) + OutputVTKHDFPoints = h5open("$(particle_savepath)_PointMeasures.vtkhdf", "w") + point_root = HDF5.create_group(OutputVTKHDFPoints, "VTKHDF") + point_positions = [to_point3d(measure.Position, T) for measure in SimMetaData.PointMeasures] + point_output_init = Vector{Any}(undef, length(point_measure_vars)) + for (index, name) in pairs(point_measure_vars) + if name in vector_fields + point_output_init[index] = fill(zero(SVector{3, T}), length(SimMetaData.PointMeasures)) + else + prop = field_map[name] + src = getproperty(SimParticles, prop) + point_output_init[index] = fill(zero(eltype(src)), length(SimMetaData.PointMeasures)) + end + end + GenerateGeometryStructure(point_root, point_measure_vars, point_output_init...; chunk_size = 1024) + GenerateStepStructure(point_root, point_measure_vars, point_output_init...) + point_handle = (file = OutputVTKHDFPoints, root = point_root) + end + + (particle_files = OutputVTKHDF, grid_files = nothing, point_files = point_handle) end end - T = eltype(eltype(SimParticles.Position)) function allocate_particle_snapshot() n = length(SimParticles.Position) positions = Vector{SVector{3, T}}(undef, n) @@ -744,8 +715,19 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, output_data[i] = similar(src, n) end end - field_data = Any[] - return ParticleSnapshot(positions, output_data, field_data) + point_positions = Vector{SVector{3, T}}(undef, length(SimMetaData.PointMeasures)) + point_output_data = Vector{Any}(undef, length(point_measure_vars)) + for (index, name) in pairs(point_measure_vars) + if name in vector_fields + point_output_data[index] = fill(zero(SVector{3, T}), length(SimMetaData.PointMeasures)) + else + prop = field_map[name] + src = getproperty(SimParticles, prop) + point_output_data[index] = fill(zero(eltype(src)), length(SimMetaData.PointMeasures)) + end + end + point_snapshot = PointMeasureSnapshot(point_positions, point_output_data) + return ParticleSnapshot(positions, output_data, point_snapshot) end function fill_particle_snapshot!(snapshot; neighbor_data = nothing) @@ -782,9 +764,11 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, copy!(buf, src) end end - FillPointMeasureData!( - snapshot.field_data, + FillPointMeasureSnapshot!( + snapshot.point_measure_snapshot.positions, + snapshot.point_measure_snapshot.output_data, SimMetaData.PointMeasures, + point_measure_vars, SimParticles, SimKernel, Dimensions; @@ -812,8 +796,6 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, snapshot.positions, output_vars, snapshot.output_data...; - field_data_names = point_measure_names, - field_data = snapshot.field_data, ) else AppendVTKHDFData( @@ -822,10 +804,28 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, snapshot.positions, output_vars, snapshot.output_data...; - field_data_names = point_measure_names, - field_data = snapshot.field_data, ) end + if !isempty(point_measure_vars) + if !SimMetaData.ExportSingleVTKHDF + SaveVTKHDF( + file_handles.point_files, + job.iteration, + point_filename(job.iteration), + snapshot.point_measure_snapshot.positions, + point_measure_vars, + snapshot.point_measure_snapshot.output_data..., + ) + else + AppendVTKHDFData( + file_handles.point_files.root, + job.time, + snapshot.point_measure_snapshot.positions, + point_measure_vars, + snapshot.point_measure_snapshot.output_data..., + ) + end + end put!(buffer_pool, snapshot) elseif job isa GridWriteJob if !SimMetaData.ExportSingleVTKHDF @@ -892,12 +892,20 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, for f in file_handles.particle_files isopen(f) && close(f) end + if file_handles.point_files !== nothing + for f in file_handles.point_files + isopen(f) && close(f) + end + end else # Close single-file handles isopen(file_handles.particle_files) && close(file_handles.particle_files) if file_handles.grid_files !== nothing isopen(file_handles.grid_files) && close(file_handles.grid_files) end + if file_handles.point_files !== nothing + isopen(file_handles.point_files.file) && close(file_handles.point_files.file) + end end end From 857776767c8cbcbfbac07a262d21a0c2931db63b Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Mon, 19 Jan 2026 00:07:15 +0100 Subject: [PATCH 5/9] Add multiblock VTKHDF option --- README.md | 4 +- src/OpenExternalPrograms.jl | 6 ++- src/ProduceHDFVTK.jl | 60 +++++++++++++++++--------- src/SimulationMetaDataConfiguration.jl | 1 + 4 files changed, 48 insertions(+), 23 deletions(-) diff --git a/README.md b/README.md index 4d678b3b..8b1411cd 100644 --- a/README.md +++ b/README.md @@ -94,7 +94,9 @@ To record a property at a specific point, populate `PointMeasures` in the `SimulationMetaData` configuration. Each `PointMeasure` names a location and a list of variables to interpolate within the neighbor stencil. These values are written into a dedicated `*_PointMeasures.vtkhdf` file so you can load probe -points alongside the particle data. +points alongside the particle data. Set `ExportMultiBlockVTKHDF=true` to emit a +single `*_MultiBlock.vtkhdf` file that groups particle data and point measures +as separate blocks for ParaView. ## Help diff --git a/src/OpenExternalPrograms.jl b/src/OpenExternalPrograms.jl index fdb51e4e..756ec2dc 100644 --- a/src/OpenExternalPrograms.jl +++ b/src/OpenExternalPrograms.jl @@ -70,7 +70,11 @@ function AutoOpenParaview(SimMetaData::SimulationMetaData, OutputVariableNames; if SimMetaData.ExportSingleVTKHDF ParaViewStateFileName = joinpath(SimMetaData.SaveLocation, SimMetaData.SimulationName) * "_SingleVTKHDFStateFile.py" - py_regex = "$(SimMetaData.SimulationName).vtkhdf" + if SimMetaData.ExportMultiBlockVTKHDF + py_regex = "$(SimMetaData.SimulationName)_MultiBlock.vtkhdf" + else + py_regex = "$(SimMetaData.SimulationName).vtkhdf" + end else ParaViewStateFileName = joinpath(SimMetaData.SaveLocation, SimMetaData.SimulationName) * "_StateFile.py" py_regex = "^$(SimMetaData.SimulationName)_(\\d+).vtk" #^ means to anchor the regex to the start of the string diff --git a/src/ProduceHDFVTK.jl b/src/ProduceHDFVTK.jl index bb17b1f1..66d62182 100644 --- a/src/ProduceHDFVTK.jl +++ b/src/ProduceHDFVTK.jl @@ -555,6 +555,7 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, particle_filename = (iter) -> "$(particle_savepath)_$(lpad(iter,6,"0")).vtkhdf" grid_filename = (iter) -> "$(grid_savepath)_$(lpad(iter,6,"0")).vtkhdf" point_filename = (iter) -> "$(particle_savepath)_PointMeasures_$(lpad(iter,6,"0")).vtkhdf" + multiblock_filename = "$(particle_savepath)_MultiBlock.vtkhdf" output_vars = SimMetaData.OutputVariables point_measure_vars = PointMeasureOutputVariableNames(SimMetaData.PointMeasures) @@ -582,13 +583,13 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, ) T = eltype(eltype(SimParticles.Position)) - to_point3d = if Dimensions == 2 - (point, ::Type{T}) -> SVector{3, T}(point[1], point[2], zero(T)) - else - (point, ::Type{T}) -> point - end - # Initialize storage for file handles + if SimMetaData.ExportMultiBlockVTKHDF && !SimMetaData.ExportSingleVTKHDF + error("ExportMultiBlockVTKHDF requires ExportSingleVTKHDF=true.") + end + + export_multiblock = SimMetaData.ExportMultiBlockVTKHDF && !isempty(point_measure_vars) + file_handles = if !SimMetaData.ExportSingleVTKHDF # Multi-file mode: vector for particle files n_outputs = if SimMetaData.OutputTimes isa AbstractVector @@ -603,7 +604,7 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, ) else # Single-file mode: handles for both files - OutputVTKHDF = h5open("$(particle_savepath).vtkhdf", "w") + OutputVTKHDF = h5open(export_multiblock ? multiblock_filename : "$(particle_savepath).vtkhdf", "w") root = HDF5.create_group(OutputVTKHDF, "VTKHDF") output_data_init = Vector{Any}(undef, length(output_vars)) @@ -619,14 +620,25 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, end end + root_particles = root + root_points = nothing + + if export_multiblock + HDF5.attrs(root)["Version"] = Int32.([2, 3]) + write_ascii_attribute(root, "Type", "MultiBlock") + blocks = HDF5.create_group(root, "Blocks") + root_particles = HDF5.create_group(blocks, "Particles") + root_points = HDF5.create_group(blocks, "PointMeasures") + end + GenerateGeometryStructure( - root, + root_particles, output_vars, output_data_init...; chunk_size = 1024, ) GenerateStepStructure( - root, + root_particles, output_vars, output_data_init...; ) @@ -653,9 +665,11 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, point_handle = nothing if !isempty(point_measure_vars) - OutputVTKHDFPoints = h5open("$(particle_savepath)_PointMeasures.vtkhdf", "w") - point_root = HDF5.create_group(OutputVTKHDFPoints, "VTKHDF") - point_positions = [to_point3d(measure.Position, T) for measure in SimMetaData.PointMeasures] + point_root = root_points + if !export_multiblock + OutputVTKHDFPoints = h5open("$(particle_savepath)_PointMeasures.vtkhdf", "w") + point_root = HDF5.create_group(OutputVTKHDFPoints, "VTKHDF") + end point_output_init = Vector{Any}(undef, length(point_measure_vars)) for (index, name) in pairs(point_measure_vars) if name in vector_fields @@ -668,16 +682,18 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, end GenerateGeometryStructure(point_root, point_measure_vars, point_output_init...; chunk_size = 1024) GenerateStepStructure(point_root, point_measure_vars, point_output_init...) - point_handle = (file = OutputVTKHDFPoints, root = point_root) + point_handle = export_multiblock ? (file = OutputVTKHDF, root = point_root, shared = true) : (file = OutputVTKHDFPoints, root = point_root, shared = false) end - (particle_files = OutputVTKHDF, grid_files = OutputVTKHDFGrid, point_files = point_handle) + (particle_files = OutputVTKHDF, grid_files = OutputVTKHDFGrid, point_files = point_handle, particle_root = root_particles) else point_handle = nothing if !isempty(point_measure_vars) - OutputVTKHDFPoints = h5open("$(particle_savepath)_PointMeasures.vtkhdf", "w") - point_root = HDF5.create_group(OutputVTKHDFPoints, "VTKHDF") - point_positions = [to_point3d(measure.Position, T) for measure in SimMetaData.PointMeasures] + point_root = root_points + if !export_multiblock + OutputVTKHDFPoints = h5open("$(particle_savepath)_PointMeasures.vtkhdf", "w") + point_root = HDF5.create_group(OutputVTKHDFPoints, "VTKHDF") + end point_output_init = Vector{Any}(undef, length(point_measure_vars)) for (index, name) in pairs(point_measure_vars) if name in vector_fields @@ -690,10 +706,10 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, end GenerateGeometryStructure(point_root, point_measure_vars, point_output_init...; chunk_size = 1024) GenerateStepStructure(point_root, point_measure_vars, point_output_init...) - point_handle = (file = OutputVTKHDFPoints, root = point_root) + point_handle = export_multiblock ? (file = OutputVTKHDF, root = point_root, shared = true) : (file = OutputVTKHDFPoints, root = point_root, shared = false) end - (particle_files = OutputVTKHDF, grid_files = nothing, point_files = point_handle) + (particle_files = OutputVTKHDF, grid_files = nothing, point_files = point_handle, particle_root = root_particles) end end @@ -799,7 +815,7 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, ) else AppendVTKHDFData( - root, + file_handles.particle_root, job.time, snapshot.positions, output_vars, @@ -904,7 +920,9 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, isopen(file_handles.grid_files) && close(file_handles.grid_files) end if file_handles.point_files !== nothing - isopen(file_handles.point_files.file) && close(file_handles.point_files.file) + if !file_handles.point_files.shared + isopen(file_handles.point_files.file) && close(file_handles.point_files.file) + end end end end diff --git a/src/SimulationMetaDataConfiguration.jl b/src/SimulationMetaDataConfiguration.jl index 437950b2..ba107d55 100644 --- a/src/SimulationMetaDataConfiguration.jl +++ b/src/SimulationMetaDataConfiguration.jl @@ -48,6 +48,7 @@ struct StoreLog <: LogMode end IndexCounter::Int = 0 VisualizeInParaview::Bool = true ExportSingleVTKHDF::Bool = true + ExportMultiBlockVTKHDF::Bool = false ExportGridCells::Bool = false ExportGridCellParticleCounts::Bool = false OutputVariables::Vector{String} = [ From 2a08b2219a954cf7ad160b24841da76ab87d9604 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Mon, 19 Jan 2026 00:17:36 +0100 Subject: [PATCH 6/9] Handle empty probe neighborhoods --- src/PointMeasurements.jl | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/src/PointMeasurements.jl b/src/PointMeasurements.jl index ee06f1f9..8cdfb227 100644 --- a/src/PointMeasurements.jl +++ b/src/PointMeasurements.jl @@ -152,7 +152,15 @@ function AppendPointMeasureResults!(FieldData, Measure::PointMeasure{D, T}, end if !HasCandidates - error("PointMeasure found no neighbor candidates; adjust probe position or cell list.") + for (Index, Variable) in pairs(Variables) + if Variable in vector_fields + push!(FieldData, MissingVectorValue(Val(Dimensions), T)) + else + Source = Sources[Index] + push!(FieldData, MissingScalarValue(eltype(Source))) + end + end + return FieldData end for (Index, Variable) in pairs(Variables) @@ -173,6 +181,18 @@ function AppendPointMeasureResults!(FieldData, Measure::PointMeasure{D, T}, return FieldData end +@inline function MissingScalarValue(::Type{T}) where {T} + return T <: AbstractFloat ? T(NaN) : zero(T) +end + +@inline function MissingVectorValue(::Val{2}, ::Type{T}) where {T} + return SVector{3, T}(MissingScalarValue(T), MissingScalarValue(T), MissingScalarValue(T)) +end + +@inline function MissingVectorValue(::Val{3}, ::Type{T}) where {T} + return SVector{3, T}(MissingScalarValue(T), MissingScalarValue(T), MissingScalarValue(T)) +end + function FillPointMeasureSnapshot!(snapshot_positions, snapshot_output_data, PointMeasures::Vector{<:PointMeasure}, OutputVariables::Vector{String}, From e6301bd48cbb2aafc57442858bf42d767cb14d2e Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Mon, 19 Jan 2026 00:36:06 +0100 Subject: [PATCH 7/9] Allow multiblock without probes --- src/ProduceHDFVTK.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/ProduceHDFVTK.jl b/src/ProduceHDFVTK.jl index 66d62182..04183c9b 100644 --- a/src/ProduceHDFVTK.jl +++ b/src/ProduceHDFVTK.jl @@ -588,7 +588,7 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, error("ExportMultiBlockVTKHDF requires ExportSingleVTKHDF=true.") end - export_multiblock = SimMetaData.ExportMultiBlockVTKHDF && !isempty(point_measure_vars) + export_multiblock = SimMetaData.ExportMultiBlockVTKHDF file_handles = if !SimMetaData.ExportSingleVTKHDF # Multi-file mode: vector for particle files @@ -605,8 +605,12 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, else # Single-file mode: handles for both files OutputVTKHDF = h5open(export_multiblock ? multiblock_filename : "$(particle_savepath).vtkhdf", "w") - root = HDF5.create_group(OutputVTKHDF, "VTKHDF") - + if !isempty(point_measure_vars) + root_points = create_ordered_group(root, "PointMeasures") + end + if !isempty(point_measure_vars) + create_soft_link("/VTKHDF/PointMeasures", assembly, "PointMeasures") + end output_data_init = Vector{Any}(undef, length(output_vars)) for (i, name) in pairs(output_vars) prop = get(field_map, name, nothing) From ca31137bead752e0696afdeae46b743c3d3b27be Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Mon, 19 Jan 2026 00:44:26 +0100 Subject: [PATCH 8/9] Inline multiblock helpers --- src/ProduceHDFVTK.jl | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/ProduceHDFVTK.jl b/src/ProduceHDFVTK.jl index 04183c9b..c3698378 100644 --- a/src/ProduceHDFVTK.jl +++ b/src/ProduceHDFVTK.jl @@ -63,21 +63,6 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, HDF5.write_attribute(attr, dtype, value) end - """ - Compute points and connectivity information for an unstructured grid given - `UniqueCells` from the SPH cell list. Returns `(points, connectivity, - offsets, cell_types, cell_ids, dims)` where `dims` is 2 or 3. - """ - function compute_grid_geometry(SimKernel, UniqueCells) - ExtractDimensionality(::AbstractVector{CartesianIndex{N}}) where N = N - - dims = ExtractDimensionality(UniqueCells) - - dx = dy = dz = SimKernel.H - - if dims == 2 - minx, maxx = minimum(ci -> ci[1], UniqueCells), maximum(ci -> ci[1], UniqueCells) - miny, maxy = minimum(ci -> ci[2], UniqueCells), maximum(ci -> ci[2], UniqueCells) nx = maxx - minx + 1 elseif dims == 3 minx, maxx = minimum(ci -> ci[1], UniqueCells), maximum(ci -> ci[1], UniqueCells) @@ -556,6 +541,21 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, grid_filename = (iter) -> "$(grid_savepath)_$(lpad(iter,6,"0")).vtkhdf" point_filename = (iter) -> "$(particle_savepath)_PointMeasures_$(lpad(iter,6,"0")).vtkhdf" multiblock_filename = "$(particle_savepath)_MultiBlock.vtkhdf" + + function ordered_group(parent, name::String) + gcpl = HDF5.API.h5p_create(HDF5.API.H5P_GROUP_CREATE) + flags = HDF5.API.H5P_CRT_ORDER_TRACKED | HDF5.API.H5P_CRT_ORDER_INDEXED + HDF5.API.h5pset_link_creation_order(gcpl, flags) + HDF5.API.h5pset_attr_creation_order(gcpl, flags) + group_id = HDF5.API.h5gcreate(parent.id, name, HDF5.API.H5P_DEFAULT, gcpl, HDF5.API.H5P_DEFAULT) + HDF5.API.h5p_close(gcpl) + return HDF5.Group(group_id) + end + + function create_soft_link(target_path::String, parent, name::String) + HDF5.API.h5l_create_soft(target_path, parent.id, name, HDF5.API.H5P_DEFAULT, HDF5.API.H5P_DEFAULT) + return nothing + end output_vars = SimMetaData.OutputVariables point_measure_vars = PointMeasureOutputVariableNames(SimMetaData.PointMeasures) @@ -605,10 +605,10 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, else # Single-file mode: handles for both files OutputVTKHDF = h5open(export_multiblock ? multiblock_filename : "$(particle_savepath).vtkhdf", "w") - if !isempty(point_measure_vars) - root_points = create_ordered_group(root, "PointMeasures") - end - if !isempty(point_measure_vars) + root = export_multiblock ? ordered_group(OutputVTKHDF, "VTKHDF") : HDF5.create_group(OutputVTKHDF, "VTKHDF") + root_particles = ordered_group(root, "Particles") + root_points = ordered_group(root, "PointMeasures") + assembly = ordered_group(root, "Assembly") create_soft_link("/VTKHDF/PointMeasures", assembly, "PointMeasures") end output_data_init = Vector{Any}(undef, length(output_vars)) From aa6031a60d33dfa23ba6a300dec303f7a3b0ac46 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Mon, 19 Jan 2026 00:48:44 +0100 Subject: [PATCH 9/9] Restore multiblock helpers --- src/ProduceHDFVTK.jl | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/ProduceHDFVTK.jl b/src/ProduceHDFVTK.jl index c3698378..de9e566b 100644 --- a/src/ProduceHDFVTK.jl +++ b/src/ProduceHDFVTK.jl @@ -63,6 +63,21 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, HDF5.write_attribute(attr, dtype, value) end + function create_ordered_group(parent, name::String) + gcpl = HDF5.API.h5p_create(HDF5.API.H5P_GROUP_CREATE) + flags = HDF5.API.H5P_CRT_ORDER_TRACKED | HDF5.API.H5P_CRT_ORDER_INDEXED + HDF5.API.h5pset_link_creation_order(gcpl, flags) + HDF5.API.h5pset_attr_creation_order(gcpl, flags) + group_id = HDF5.API.h5gcreate(parent.id, name, HDF5.API.H5P_DEFAULT, gcpl, HDF5.API.H5P_DEFAULT) + HDF5.API.h5p_close(gcpl) + return HDF5.Group(group_id) + end + + function create_soft_link(target_path::String, parent, name::String) + HDF5.API.h5l_create_soft(target_path, parent.id, name, HDF5.API.H5P_DEFAULT, HDF5.API.H5P_DEFAULT) + return nothing + end + nx = maxx - minx + 1 elseif dims == 3 minx, maxx = minimum(ci -> ci[1], UniqueCells), maximum(ci -> ci[1], UniqueCells) @@ -542,21 +557,6 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, point_filename = (iter) -> "$(particle_savepath)_PointMeasures_$(lpad(iter,6,"0")).vtkhdf" multiblock_filename = "$(particle_savepath)_MultiBlock.vtkhdf" - function ordered_group(parent, name::String) - gcpl = HDF5.API.h5p_create(HDF5.API.H5P_GROUP_CREATE) - flags = HDF5.API.H5P_CRT_ORDER_TRACKED | HDF5.API.H5P_CRT_ORDER_INDEXED - HDF5.API.h5pset_link_creation_order(gcpl, flags) - HDF5.API.h5pset_attr_creation_order(gcpl, flags) - group_id = HDF5.API.h5gcreate(parent.id, name, HDF5.API.H5P_DEFAULT, gcpl, HDF5.API.H5P_DEFAULT) - HDF5.API.h5p_close(gcpl) - return HDF5.Group(group_id) - end - - function create_soft_link(target_path::String, parent, name::String) - HDF5.API.h5l_create_soft(target_path, parent.id, name, HDF5.API.H5P_DEFAULT, HDF5.API.H5P_DEFAULT) - return nothing - end - output_vars = SimMetaData.OutputVariables point_measure_vars = PointMeasureOutputVariableNames(SimMetaData.PointMeasures) @@ -605,10 +605,10 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, else # Single-file mode: handles for both files OutputVTKHDF = h5open(export_multiblock ? multiblock_filename : "$(particle_savepath).vtkhdf", "w") - root = export_multiblock ? ordered_group(OutputVTKHDF, "VTKHDF") : HDF5.create_group(OutputVTKHDF, "VTKHDF") - root_particles = ordered_group(root, "Particles") - root_points = ordered_group(root, "PointMeasures") - assembly = ordered_group(root, "Assembly") + root = export_multiblock ? create_ordered_group(OutputVTKHDF, "VTKHDF") : HDF5.create_group(OutputVTKHDF, "VTKHDF") + root_particles = create_ordered_group(root, "Particles") + root_points = create_ordered_group(root, "PointMeasures") + assembly = create_ordered_group(root, "Assembly") create_soft_link("/VTKHDF/PointMeasures", assembly, "PointMeasures") end output_data_init = Vector{Any}(undef, length(output_vars))