Read VTK

Project.toml

@@ -22,6 +22,7 @@ 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"

examples/readvtk/readvtk.jl

@@ -0,0 +1,18 @@
+# Example File to read a vtk file and convert it to InitialCondition
+using TrixiParticles
+
+# Create a example rectangle
+rectangle = RectangularShape(0.1, (10, 10), (0, 0), density=1.5, velocity=(1.0, -2.0),
+                             pressure=1000.0)
+
+trixi2vtk(rectangle; filename="rectangle", output_directory="out",
+          custom_quantity=nothing)
+
+filename = "rectangle"
+file = joinpath("out", filename * ".vtu")
+
+# Read the vtk file and convert it to an 'InitialCondition'
+ic = vtk2trixi(file)
+
+trixi2vtk(ic; filename="readvtk", output_directory="out",
+          custom_quantity=nothing)

src/TrixiParticles.jl

@@ -19,6 +19,7 @@ using LinearAlgebra: norm, dot, I, tr, inv, pinv, det
 using MuladdMacro: @muladd
 using Polyester: Polyester, @batch
 using Printf: @printf, @sprintf
+using ReadVTK: ReadVTK
 using RecipesBase: RecipesBase, @series
 using Random: seed!
 using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!,
@@ -75,7 +76,7 @@ export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureEx
 
 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,49 @@
+"""
+    vtk2trixi(file::String)
+
+Convert data from VTK-file to InitialCondition 
+
+# Arguments
+- `file`: Name of the file to be loaded. 
+"""
+function vtk2trixi(file)
+    vtk_file = ReadVTK.VTKFile(file)
+
+    # Retrieve data fields (e.g., pressure, velocity, ...)
+    point_data = ReadVTK.get_point_data(vtk_file)
+    field_data = ReadVTK.get_field_data(vtk_file)
+
+    # Retrieve fields
+    ndims = ReadVTK.get_data(field_data["ndims"])
+    coordinates = ReadVTK.get_points(vtk_file)[1:ndims[1], :]
+
+    fields = ["velocity", "density", "pressure", "mass"]
+    results = Dict{String, Array{Float64}}()
+
+    for field in fields
+        found = false
+        for k in ReadVTK.keys(point_data)
+            if !isnothing(match(Regex("$field"), k))
+                results[field] = ReadVTK.get_data(point_data[k])
+                found = true
+                break
+            end
+        end
+        if !found
+            # Set fields to zero if not found
+            if field in ["density", "pressure", "mass"]
+                results[field] = zeros(size(coordinates, 2))
+            else
+                results[field] = zero(coordinates)
+            end
+            @info "No '$field' field found in VTK file. Will be set to zero."
+        end
+    end
+
+    return InitialCondition(
+                            ; coordinates,
+                            velocity=results["velocity"],
+                            mass=results["mass"],
+                            density=results["density"],
+                            pressure=results["pressure"])
+end

src/visualization/write2vtk.jl

@@ -123,6 +123,7 @@ function trixi2vtk(v_, u_, t, system_, periodic_box; output_directory="out", pre
         # Store particle index
         vtk["index"] = active_particles(system)
         vtk["time"] = t
+        vtk["ndims"] = ndims(system)
 
         if write_meta_data
             vtk["solver_version"] = git_hash
@@ -194,6 +195,7 @@ function trixi2vtk(coordinates; output_directory="out", prefix="", filename="coo
     vtk_grid(file, points, cells) do vtk
         # Store particle index.
         vtk["index"] = [i for i in axes(coordinates, 2)]
+        vtk["ndims"] = size(coordinates, 1)
 
         # Extract custom quantities for this system.
         for (key, quantity) in custom_quantities
@@ -202,7 +204,6 @@ function trixi2vtk(coordinates; output_directory="out", prefix="", filename="coo
             end
         end
     end
-
     return file
 end
 

test/preprocessing/preprocessing.jl

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

test/preprocessing/readvtk/vtk2trixi.jl

@@ -0,0 +1,40 @@
+@testset verbose=true "vtk2trixi" begin
+    @testset verbose=true "Functionality Check - 2D" begin
+
+        # 'InitialCondition'-Files
+        saved_ic = RectangularShape(0.1, (10, 20),
+                                    (0, 0), density=1.5,
+                                    velocity=(1.0, 2.0), pressure=1000.0)
+
+        file = trixi2vtk(saved_ic; output_directory="test/preprocessing/readvtk",
+                         filename="initial_condition")
+
+        loaded_ic = vtk2trixi(joinpath("test/preprocessing/readvtk",
+                                       "initial_condition.vtu"))
+
+        @test isapprox(saved_ic.coordinates, loaded_ic.coordinates)
+        @test isapprox(saved_ic.velocity, loaded_ic.velocity)
+        @test isapprox(saved_ic.density, loaded_ic.density)
+        @test isapprox(saved_ic.pressure, loaded_ic.pressure)
+        @test isapprox(saved_ic.mass, loaded_ic.mass)
+
+        # 'Fluidsystem'-Files
+        state_equation = StateEquationCole(; sound_speed=1.0, reference_density=1.0,
+                                           exponent=1, clip_negative_pressure=false)
+        fluid_system = WeaklyCompressibleSPHSystem(saved_ic,
+                                                   SummationDensity(), state_equation,
+                                                   WendlandC2Kernel{2}(), 0.2)
+
+        saved_fs = trixi2vtk(saved_ic.velocity, saved_ic.coordinates, 0.0, fluid_system,
+                             nothing;
+                             output_directory="test/preprocessing/readvtk", iter=1)
+
+        loaded_fs = vtk2trixi(joinpath("test/preprocessing/readvtk", "fluid_1.vtu"))
+
+        @test isapprox(saved_fs.coordinates, loaded_fs.coordinates)
+        @test isapprox(saved_fs.velocity, loaded_fs.velocity)
+        @test isapprox(saved_fs.density, loaded_fs.density)
+        @test isapprox(saved_fs.pressure, loaded_fs.pressure)
+        @test isapprox(saved_fs.mass, loaded_fs.mass)
+    end
+end