add vtk routine to save marker chain

.github/workflows/Documenter.yml

@@ -33,7 +33,7 @@ jobs:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
         with:
-          version: '1.10'
+          version: '1.11'
       - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@latest
       - name: Install dependencies

Project.toml

@@ -38,7 +38,7 @@ GeoParams = "0.6.7"
 HDF5 = "0.17.1"
 ImplicitGlobalGrid = "0.15"
 JLD2 = "0.4, 0.5"
-JustPIC = "0.5.3"
+JustPIC = "0.5"
 MPI = "0.20"
 MuladdMacro = "0.2"
 ParallelStencil = "0.13.6, 0.14"

src/IO/DataIO.jl

@@ -4,6 +4,7 @@ using WriteVTK
 using HDF5
 using JLD2
 using MPI
+using StaticArrays
 
 import ..JustRelax: Geometry
 import ..JustRelax: IGG
@@ -24,7 +25,7 @@ export checkpointing_jld2, load_checkpoint_jld2
 
 include("VTK.jl")
 
-export VTKDataSeries, append!, save_vtk
+export VTKDataSeries, append!, save_vtk, save_marker_chain
 
 export metadata
 

src/IO/VTK.jl

@@ -122,3 +122,33 @@ function save_vtk(fname::String, xi, data::NamedTuple)
 
     return nothing
 end
+
+"""
+    save_marker_chain(fname::String, cell_vertices::LinRange{Float64}, h_vertices::Vector{Float64})
+
+Save a vector of points as a line in a VTK file. The X and Y coordinates of the points are given by `cell_vertices` and `h_vertices`, respectively.
+The Z coordinate is set to zero as we are working in 2D.
+
+## Arguments
+- `fname::String`: The name of the VTK file to save. The extension `.vtk` will be appended to the name.
+- `cell_vertices::LinRange{Float64}`: A range of X coordinates for the points.
+- `h_vertices::Vector{Float64}`: A vector of Y coordinates for the points.
+
+## Example
+```julia
+cell_vertices = LinRange(0.0, 1.0, 10)
+h_vertices = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
+save_marker_chain("Example", cell_vertices, h_vertices)
+```
+"""
+function save_marker_chain(fname::String, cell_vertices::LinRange{Float64}, h_vertices::Vector{Float64})
+    cell_vertices_vec = collect(cell_vertices)  # Convert LinRange to Vector
+    n_points = length(cell_vertices_vec)
+    points = [SVector{3, Float64}(cell_vertices_vec[i], h_vertices[i], 0.0) for i in 1:n_points]
+    lines = [MeshCell(PolyData.Lines(), 1:(n_points))]  # Create a single line connecting all points
+
+    vtk_grid(fname, points, lines) do vtk
+        vtk["Points"] = points
+    end
+    return nothing
+end

test/Project.toml

@@ -22,7 +22,7 @@ AMDGPU = "0.8, 0.9, 1"
 CUDA = "5"
 CellArrays = "0.2, 0.3"
 GeoParams = "0.6.7"
-JustPIC = "0.5.3"
+JustPIC = "0.5"
 MPI = "0.20"
 ParallelStencil = "0.11, 0.12, 0.13, 0.14"
 SpecialFunctions = "2"

test/runtests.jl

@@ -55,15 +55,15 @@ function runtests()
             try
                 @testset "$(basename(f))" begin
                     n = 2
-                    p= run(`$(mpiexec()) -n $n $(Base.julia_cmd()) -O3 --startup-file=no --check-bounds=no $(joinpath(testdir, f))`)
+                    p= run(`$(mpiexec()) -n $n $(Base.julia_cmd()) -O3 --startup-file=no $(joinpath(testdir, f))`)
                     @test success(p)
                 end
             catch ex
                 nfail += 1
             end
         else
             try
-                run(`$(Base.julia_cmd()) -O3 --startup-file=no --check-bounds=no $(joinpath(testdir, f))`)
+                run(`$(Base.julia_cmd()) -O3 --startup-file=no $(joinpath(testdir, f))`)
             catch ex
                 nfail += 1
             end

test/test_IO.jl

@@ -154,6 +154,14 @@ using WriteVTK
         DataIO.append!(vtk, (Vy=stokes.V.Vy, η=stokes.viscosity.η), dt, time)
         @test isfile(joinpath(dst, "vtk_series.pvd"))
 
+        ## Test save_marker_chain
+        nxcell, max_xcell, min_xcell = 100, 150, 75
+        initial_elevation            = 0e0
+        chain                        = init_markerchain(backend, nxcell, min_xcell, max_xcell, xvi[1], initial_elevation);
+
+        save_marker_chain(joinpath(dst,"MarkerChain"), chain.cell_vertices, chain.h_vertices)
+        @test isfile(joinpath(dst, "MarkerChain.vtp"))
+
         # 3D case
         ni = nx, ny, nz
         stokes  = StokesArrays(backend_JR, ni)

test/test_shearband2D_softening.jl

@@ -200,7 +200,7 @@ end
     @suppress begin
         iters, τII, sol = ShearBand2D()
         @test passed = iters.err_evo1[end] < 1e-6
-        @test τII[end] ≈ 1.59639 atol = 1e-4
+        @test τII[end] ≈ 1.59073 atol = 1e-4
         @test sol[end] ≈ 1.94255 atol = 1e-4
     end
 end