Read vtk

Project.toml

@@ -5,11 +5,13 @@ version = "0.2.4-dev"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FastPow = "c0e83750-1142-43a8-81cf-6c956b72b4d1"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -18,17 +20,20 @@ JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PointNeighbors = "1c4d5385-0a27-49de-8e2c-43b175c8985c"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StrideArrays = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
+VTKDataIO = "c6703add-1d23-52c6-9943-3ad88652b9b2"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
@@ -44,6 +49,7 @@ GPUArraysCore = "0.1, 0.2"
 JSON = "0.21"
 KernelAbstractions = "0.9"
 MuladdMacro = "0.2"
+OrdinaryDiffEq = "6.87.0"
 PointNeighbors = "0.4.2"
 Polyester = "0.7.10"
 RecipesBase = "1"

examples/readvtk/trixi2vtk.jl

@@ -0,0 +1,8 @@
+# Example File to read a vtk file and convert it to InitialCondition
+using TrixiParticles
+using OrdinaryDiffEq
+
+ic = vtk2trixi("examples/readvtk/fluid_1_1.vtu")
+
+trixi2vtk(ic; filename="trixi2vtk_test", output_directory="out_vtk",
+          custom_quantity=nothing)

examples/vtk-reader/particle_saving.jl

@@ -0,0 +1,61 @@
+using ReadVTK
+
+struct Particle
+    index::Int64
+    coordinates::Tuple{Float64, Float64, Float64}  # Assuming 3D coordinates
+    density::Float64
+    velocity::Tuple{Float64, Float64, Float64}  # Assuming 3D velocity
+    pressure::Float64
+end
+
+vtk_file = VTKFile("out/fluid_1_1.vtu")
+#vtk_file = VTKFile("out_vtk/rectangle_of_water.vtu")
+
+# Retrieve particle coordinates
+coords = get_points(vtk_file)
+
+# Retrieve point data fields (e.g., pressure, velocity, ...)
+vtk_point_data = get_point_data(vtk_file)
+
+# Dynamically get all available field names from point_data
+field_data = Dict()
+for field_name in keys(vtk_point_data)
+    field_data[field_name] = get_data(vtk_point_data[field_name])
+end
+
+# Create an array of Particle instances
+particles = Vector{Particle}(undef, size(coords, 2))
+
+for i in 1:size(coords, 2)
+    # Retrieve required field "index"
+    index = field_data["index"][i]
+
+    # Coordinates
+    coordinates = (coords[1, i], coords[2, i], coords[3, i])
+
+    # Retrieve each required field directly, assuming all are present
+    density = field_data["density"][i]
+    pressure = field_data["pressure"][i]
+
+    velocity = if size(field_data["velocity"], 1) == 2
+        # If velocity is 2D, add 0.0 for the z component
+        (field_data["velocity"][1, i], field_data["velocity"][2, i], 0.0)
+    else
+        # If velocity is 3D, use all three components
+        (field_data["velocity"][1, i], field_data["velocity"][2, i],
+         field_data["velocity"][3, i])
+    end
+
+    # Create a new Particle instance
+    particles[i] = Particle(index, coordinates, density, velocity, pressure)
+end
+
+# Display some particles for verification
+for particle in particles[1:2]
+    println(particle)
+end
+
+println("Coords: ", size(coords))
+println("velocity ", size(field_data["velocity"]))
+
+# TODO FieldData 

examples/vtk-reader/read_two_points_dynamic.jl

@@ -0,0 +1,23 @@
+using ReadVTK
+
+#vtk = VTKFile("out/fluid_1_0.vtu")
+#vtk = VTKFile("out_vtk/two_points.vtu")
+vtk = VTKFile("out_vtk/rectangle_of_water.vtu")
+
+point_data_file = get_point_data(vtk)
+
+point_names = keys(point_data_file)
+
+for point_name in point_names
+    point_array = point_data_file[point_name]
+
+    point_data = get_data(point_array)
+
+    println(point_name)
+    println(point_data)
+end
+
+coords = get_points(vtk)
+for i in 1:size(coords, 2)
+    println("Coords $i: (", coords[1, i], ", ", coords[2, i], ", ", coords[3, i], ")")
+end

examples/vtk-reader/read_two_points_static.jl

@@ -0,0 +1,51 @@
+using ReadVTK
+
+vtk = VTKFile("out_vtk/two_points.vtu")
+#vtk = VTKFile("out/fluid_1_0.vtu")
+
+"""
+# Fields of vtk
+
+- `filename`: original path to the VTK file that has been read in
+- `xml_file`: object that represents the XML file
+- `file_type`: currently only `"UnstructuredGrid"` or `"ImageData"` are supported
+- `version`: VTK XML file format version
+- `byte_order`: can be `LittleEndian` or `BigEndian` and must currently be the same as the system's
+- `compressor`: can be empty (no compression) or `vtkZLibDataCompressor`
+- `appended_data`: in case of appended data (see XML documentation), the data is stored here for
+                   convenient retrieval (otherwise it is empty)
+- `n_points`: number of points in the VTK file
+- `n_cells`: number of cells in the VTK file`
+"""
+
+""" Examples
+
+#cell_data = get_cell_data(vtk)
+
+point_data = get_point_data(vtk)
+# create an vector with the names of the point data like "velocity"
+names = point_data.names
+
+# create an vector with the information about the Data like "offset"
+data_arrays = point_data.data_arrays
+
+# information about the file like "UnstructedGrid"
+vtk_file = point_data.vtk_file
+
+# gets the data of the vtk-file
+appended_data = vtk.appended_data
+
+# gets the coords of all the points
+#coords = get_points(vtk)
+
+# gets the information about the cells
+#cells = get_cells(vtk)
+"""
+point_data = get_point_data(vtk)
+
+velocity_array = point_data["velocity"]
+
+velocity_data = get_data(velocity_array)
+
+println("Velocity Data:")
+println(velocity_data)

examples/vtk-reader/write_n_points.jl

@@ -0,0 +1,51 @@
+using WriteVTK
+
+# Function to create n points and their associated data
+function create_vtk_file(n::Int, coords::Matrix{Float64}, scalars::Vector{Float64},
+                         velocities::Matrix{Float64}, filename::String)
+    # Ensure the coordinates matrix has the correct size (3 x n)
+    @assert size(coords)==(3, n) "coords should be a 3 x n matrix, where n is the number of points"
+    # Ensure scalars is a vector of length n
+    @assert length(scalars)==n "scalars should be a vector of length n"
+    # Ensure velocities is a matrix of size (3 x n)
+    @assert size(velocities)==(3, n) "velocities should be a 3 x n matrix, where n is the number of points"
+
+    output_directory = "out_vtk"
+    mkpath(output_directory)
+
+    file = joinpath(output_directory, filename)
+
+    # Create VTK_VERTEX cells for all points
+    cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (i,)) for i in 1:n]
+
+    # Initialize the VTK grid with the points and cells
+    vtk_file = vtk_grid(file, coords, cells)
+
+    # Assign scalar data to the points
+    vtk_file["scalar"] = scalars
+
+    # Flatten velocities into a 1D array for VTK
+    vtk_file["velocity"] = vec(velocities)
+
+    # Save the VTK file
+    vtk_save(vtk_file)
+end
+
+# Example usage with n points
+n = 3  # Number of points
+
+# Define coordinates for n points (3 rows for x, y, z; n columns for the points)
+coords = [1.0 2.0 3.0;  # x-coordinates
+          1.0 2.0 3.0;  # y-coordinates
+          0.0 0.0 0.0]  # z-coordinates (all on the z=0 plane)
+
+# Define scalar values for each point
+scalars = [100.0, 200.0, 300.0]
+
+# Define velocity vectors for each point (3 rows for vx, vy, vz; n columns for the points)
+velocities = [10.0 0.0 5.0;  # x-components
+              0.0 5.0 2.0;  # y-components
+              0.0 0.0 0.0]  # z-components (all on the z=0 plane)
+
+# Create the VTK file
+create_vtk_file(n, coords, scalars, velocities, "n_points")

examples/vtk-reader/write_one_point.jl

@@ -0,0 +1,28 @@
+using WriteVTK
+
+output_directory = "out_vtk"
+mkpath(output_directory)
+
+file = joinpath(output_directory, "one_point")
+
+# Define the points array in the expected format: 3 rows (x, y, z), 1 column (for 1 point)
+points = [1.0  # x-coordinate
+          2.0  # y-coordinate
+          0.0] # z-coordinate
+
+# Create a VTK_VERTEX cell for the single point
+cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (1,))]  # Cell for the single point
+
+# Initialize the VTK grid with the point and cell
+vtk_file = vtk_grid(file, points, cells)
+
+# Assign scalar data to the point
+vtk_file["scalar"] = [100.0]  # Scalar value for the point
+
+# Assign vector data (e.g., velocity) to the point
+vtk_file["velocity"] = [10.0, 0.0, 0.0]  # Velocity for the point
+
+# Save the VTK file
+vtk_save(vtk_file)
+
+# Surprisingly, the code works not for one point, but for >=2 points.

examples/vtk-reader/write_two_points.jl

@@ -0,0 +1,28 @@
+using WriteVTK
+
+output_directory = "out_vtk"
+mkpath(output_directory)
+
+file = joinpath(output_directory, "two_points")
+
+# Define the points array in the expected format: 3 rows (x, y, z), 2 columns (for 2 points)
+points = [1.0 3.0  # x-coordinates
+          2.0 4.0  # y-coordinates
+          0.0 0.0] # z-coordinates (both points on the z=0 plane)
+
+# Create VTK_VERTEX cells for both points
+cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (1,)),  # Cell for first point
+    MeshCell(VTKCellTypes.VTK_VERTEX, (2,))]  # Cell for second point
+
+# Initialize the VTK grid with the points and cells
+vtk_file = vtk_grid(file, points, cells)
+
+# Assign scalar data to the two points
+vtk_file["scalar"] = [100.0, 200.0]  # Scalar values for the two points
+
+# Assign vector data (e.g., velocity) to the two points
+vtk_file["velocity"] = [1.0, 2.0, 3.0,   # Velocity for first point
+    4.0, 5.0, 6.0]    # Velocity for second point
+
+# Save the VTK file
+vtk_save(vtk_file)

examples/vtk-reader/write_with_trixi.jl

@@ -0,0 +1,8 @@
+using TrixiParticles
+
+rectangle = RectangularShape(0.1, (1, 2), (0.0, 0.0), velocity=[10.0, 0.0], density=1000.0)
+
+custom_quantity = ones(nparticles(rectangle))
+
+trixi2vtk(rectangle; filename="rectangle_of_water", output_directory="out_vtk",
+          custom_quantity=nothing)

src/TrixiParticles.jl

@@ -72,7 +72,7 @@ export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureEx
        PressureMirroring, PressureZeroing, BoundaryModelLastiwka
 export BoundaryMovement
 export examples_dir, validation_dir, trixi_include
-export trixi2vtk
+export trixi2vtk, vtk2trixi
 export RectangularTank, RectangularShape, SphereShape, ComplexShape
 export WindingNumberHormann, WindingNumberJacobson
 export VoxelSphere, RoundSphere, reset_wall!, extrude_geometry, load_geometry

src/preprocessing/preprocessing.jl

@@ -1,2 +1,3 @@
 include("geometries/geometries.jl")
 include("point_in_poly/point_in_poly.jl")
+include("readvtk/vtk2trixi.jl")

src/preprocessing/readvtk/vtk2trixi.jl

@@ -0,0 +1,36 @@
+using ReadVTK
+
+"""
+    Convert data from VTK-file to InitialCondition
+
+    # FluidSystem data only
+"""
+
+function vtk2trixi(filename)
+    vtk_file = VTKFile(filename)
+
+    # Retrieve particle coordinates
+    coordinates = get_points(vtk_file)
+
+    # Retrieve point data fields (e.g., pressure, velocity, ...)
+    vtk_point_data = get_point_data(vtk_file)
+
+    # create field data arrays
+    density = get_data(vtk_point_data["density"])
+
+    pressure = get_data(vtk_point_data["pressure"])
+
+    velocity = get_data(vtk_point_data["velocity"])
+    if size(velocity, 1) == 2
+        # If velocity is 2D, add 0.0 for the z component
+        velocity = vcat(velocity, zeros(1, size(velocity, 2)))
+    end
+
+    mass = ones(size(coordinates, 2))
+
+    return InitialCondition(; coordinates, velocity, mass, density, pressure)
+end
+
+# TODO: edit the mass array --> InitialCondition needs a mass or a particle_spacing
+# TODO: example file in folder examples/readvtk
+# TODO: make it work with 2D velocity --> In ParaView the velocity vecor is 3D after trixi2vtk. It should be 2D if the initial data is 2D.