Merge pull request #83 from JuliaGeodynamics/pa-pvd

Paraview collection

Author: aelligp

docs/make.jl

@@ -108,6 +108,7 @@ makedocs(;
             "Projection" =>  "man/projection.md",
             "Surfaces" =>  "man/surfaces.md",
             "Paraview output" => "man/paraview_output.md",
+            "Paraview collection" => "man/paraview_collection.md",
             "Tools" => "man/tools.md",
             "Visualisation" => "man/visualise.md",
             "Gravity code" => "man/gravity_code.md",

docs/src/man/paraview_collection.md

@@ -0,0 +1,6 @@
+# Paraview collection
+
+We have one main routine to generate `*.pvd` files from existing `vtk` files. This is useful if no `*.pvd` file was generated during the simulation, or if you want to generate a `*.pvd` file from a collection of `vtk` files that were generated in different simulations. The `*.pvd` file can be used to animate temporal data in paraview. You can either provide a Vector of the files or the specific time step or the function reads the directory and assigns a pseudo time step to the `*.pvd` file.
+```@docs
+GeophysicalModelGenerator.make_paraview_collection
+```

src/GeophysicalModelGenerator.jl

@@ -36,6 +36,7 @@ include("coord_conversion.jl")
 include("nearest_points.jl")
 include("utils.jl")
 include("Paraview_output.jl")
+include("Paraview_collection.jl")
 include("transformation.jl")
 include("voxel_gravity.jl")
 include("LaMEM_io.jl")
@@ -52,13 +53,13 @@ include("movies_from_pics.jl")
 # GMT routines
 
 """
-        ImportTopo     
+        ImportTopo
 Optional routine that imports topography. It requires you to load `GMT`
 """
 function ImportTopo end
 
 """
-        ImportGeoTIFF     
+        ImportGeoTIFF
 Optional routine that imports GeoTIFF images. It requires you to load `GMT`
 """
 function ImportGeoTIFF end

src/Paraview_collection.jl

@@ -0,0 +1,95 @@
+using LightXML
+using WriteVTK, Printf
+
+export make_paraview_collection
+
+"""
+    make_paraview_collection(; dir=pwd(), pvd_name=nothing, files=nothing, file_extension = ".vts", time = nothing)
+
+In case one has a list of `*.vtk` files, this routine creates a `*.pvd` file that can be opened in Paraview.
+This is useful if you previously saved vtk files but didnt save it as a collection in the code itself.
+
+Optional options
+===
+- `dir`:    directory where the `*.vtk` are stored.
+- `pvd_name`:  filename of the resulting `*.pvd` file without extension; if not specified, `full_simulation` is used.
+- `files`:  Vector of the `*.vtk` files without extension; if not specified, all `*.vtk` files in the directory are used.
+- `file_extension`:  file extension of the vtk files. Default is `.vts` but all `vt*` work.
+- `time`:  Vector of the timesteps; if not specified, pseudo time steps are assigned.
+"""
+function make_paraview_collection(; dir=pwd(), pvd_name=nothing, files=nothing, file_extension = ".vts", time = nothing)
+
+    # if no files are given, use all vtm files in the directory
+    curdir = pwd()
+    cd(dir)
+    if isnothing(files)
+        files = filter(endswith(file_extension), readdir())
+        files = joinpath.(dir, files)
+    end
+    # if no time is given, use pseudo time steps
+    if isnothing(time)
+        time = [string(i) for i in 1:length(files)]
+    end
+    # Check that the arrays have the same length
+    @assert length(time) == length(files)
+    @assert length(files) > 0
+
+    # Create a new paraview collection
+    if isnothing(pvd_name)
+        pvd_name = "full_simulation"
+    end
+
+    pvd_name = joinpath(curdir, pvd_name)
+
+    make_paraview_collection(pvd_name, files, time)
+    cd(curdir)  # return to original directory
+    return pvd_name
+end
+
+"""
+    make_paraview_collection(pvd_name::String, files::Vector{String}, time::Vector{String})
+This function makes a `*.pvd` file from a provided list of `*.vtk` files and time steps.
+
+Inputs are:
+- `pvd_name`:  filename of the resulting `*.pvd` file without extension.
+- `files`: Vector of desired files to be included in the `*.pvd` file.
+- `time`: Vector of time steps corresponding to the `files` vector.
+"""
+
+function make_paraview_collection(pvd_name::String, files::Vector{String}, time::Vector{String})
+    # Check that the arrays have the same length
+    @assert length(time) == length(files)
+
+    # Create a new paraview collection
+    pvd = paraview_collection(pvd_name)
+    # save the collection
+    vtk_save(pvd)
+
+    # Parse the XML file
+    xdoc = parse_file("$pvd_name.pvd")
+
+    # Find the <Collection> element
+    collection = find_element(root(xdoc), "Collection")
+
+    # Loop over the time and data values
+    for i in 1:length(time)
+        # Create a new <DataSet> element
+        new_dataset = new_child(collection, "DataSet")
+        # Set the attributes of the new <DataSet> element
+        set_attribute(new_dataset, "timestep", time[i])
+        set_attribute(new_dataset, "part", "0")
+        set_attribute(new_dataset, "file", files[i])
+    end
+
+    # Write the modified XML to a new file of your choice (or overwrite the original)
+    open("$pvd_name.pvd", "w") do f
+        print(f, xdoc)
+    end
+    # pretty print the xml file
+    str = read("$pvd_name.pvd", String)
+    str_out = replace(str, "/>" => "/>\n")
+    open("$pvd_name.pvd", "w") do f
+        print(f, str_out)
+    end
+    return pvd_name
+end

test/runtests.jl

@@ -10,6 +10,9 @@ end
 @testset "Paraview" begin
     include("test_paraview.jl")
 end
+@testset "Paraview collection" begin
+    include("test_paraview_collection.jl")
+end
 @testset "Gravity model" begin
     include("test_voxel_gravity.jl")
 end
@@ -57,8 +60,7 @@ end
     include("test_create_movie.jl")
 end
 
-# Cleanup 
+# Cleanup
 foreach(rm, filter(endswith(".vts"), readdir()))
 foreach(rm, filter(endswith(".vtu"), readdir()))
 rm("./markers/",recursive=true)
-

test/test_paraview_collection.jl

@@ -0,0 +1,52 @@
+using Test
+using GeophysicalModelGenerator, WriteVTK
+@testset "Paraview collection" begin
+
+x, y, z = 0:10, 1:6, 2:0.1:3
+times = range(0, 1; step = 1)
+
+#generate `*.vti` files
+for (n, time) ∈ enumerate(times)
+    vtk_grid("./test_files/test_vti_$n", x, y, z) do vtk
+        vtk["Pressure"] = rand(length(x), length(y), length(z))
+    end
+end
+
+# Generate a 3D grid
+Lon,Lat,Depth   =   LonLatDepthGrid(10:20,30:40,(-300:25:0)km);
+Data            =   Depth*2; # some data
+Data_set        =   GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon))
+Write_Paraview(Data_set, "./test_files/test_depth3D")
+
+make_paraview_collection(;dir = "./test_files", pvd_name="test", file_extension=".vti")
+@test isfile("test.pvd")
+@test filesize("test.pvd") == 317
+
+make_paraview_collection(;dir = "./test_files", file_extension=".vti")
+@test isfile("full_simulation.pvd")
+@test filesize("full_simulation.pvd") == 317
+
+make_paraview_collection(;dir = "./test_files")
+@test isfile("full_simulation.pvd")
+@test filesize("full_simulation.pvd") == 251
+
+
+files = ["test_files/test_vti_1.vti", "test_files/test_vti_2.vti"]
+time  = ["1.0", "2.0"]
+make_paraview_collection("test2", files, time)
+@test isfile("test2.pvd")
+@test filesize("test2.pvd") == 317
+
+make_paraview_collection(; pvd_name="test3", files=files, time=time)
+@test isfile("test3.pvd")
+@test filesize("test3.pvd") == 317
+
+rm("test.pvd")
+rm("full_simulation.pvd")
+rm("test_files/test_depth3D.vts")
+rm("test_files/test_vti_1.vti")
+rm("test_files/test_vti_2.vti")
+rm("test2.pvd")
+rm("test3.pvd")
+
+end