Granular.jl

Julia package for granular dynamics simulation
git clone git://src.adamsgaard.dk/Granular.jl # fast
git clone https://src.adamsgaard.dk/Granular.jl.git # slow
Log | Files | Refs | README | LICENSE Back to index

io.jl (58171B)


      1 import WriteVTK
      2 import Pkg
      3 import Dates
      4 import JLD2
      5 using DelimitedFiles
      6 
      7 ## IO functions
      8 
      9 export writeSimulation
     10 """
     11     writeSimulation(simulation::Simulation;
     12                          filename::String="",
     13                          folder::String=".",
     14                          verbose::Bool=true)
     15 
     16 Write all content from `Simulation` to disk in JLD2 format.  If the `filename` 
     17 parameter is not specified, it will be saved to a subdirectory under the current 
     18 directory named after the simulation identifier `simulation.id`.
     19 """
     20 function writeSimulation(simulation::Simulation;
     21                          filename::String="",
     22                          folder::String=".",
     23                          verbose::Bool=true)
     24     if filename == ""
     25         folder = folder * "/" * simulation.id
     26         mkpath(folder)
     27         filename = string(folder, "/", simulation.id, ".",
     28                           simulation.file_number, ".jld2")
     29     end
     30 
     31     JLD2.@save(filename, simulation)
     32 
     33     if verbose
     34         @info "simulation written to $filename"
     35     end
     36     nothing
     37 end
     38 
     39 export readSimulation
     40 """
     41     readSimulation(filename::String="";
     42                    verbose::Bool=true)
     43 
     44 Return `Simulation` content read from disk using the JLD2 format.
     45 
     46 # Arguments
     47 * `filename::String`: path to file on disk containing the simulation
     48     information.
     49 * `verbose::Bool=true`: confirm to console that the file has been read.
     50 """
     51 function readSimulation(filename::String;
     52                          verbose::Bool=true)
     53     simulation = createSimulation()
     54     JLD2.@load(filename, simulation)
     55     return simulation
     56 end
     57 """
     58     readSimulation(simulation::Simulation;
     59                    step::Integer = -1,
     60                    verbose::Bool = true)
     61 
     62 Read the simulation state from disk and return as new simulation object.
     63 
     64 # Arguments
     65 * `simulation::Simulation`: use the `simulation.id` to determine the file name
     66     to read from, and read information from the file into this object.
     67 * `step::Integer=-1`: attempt to read this output file step. At its default
     68     value (`-1`), the function will try to read the latest file, determined by
     69     calling [`readSimulationStatus`](@ref).
     70 * `verbose::Bool=true`: confirm to console that the file has been read.
     71 """
     72 function readSimulation(simulation::Simulation;
     73                          step::Integer = -1,
     74                          verbose::Bool = true)
     75     if step == -1
     76         step = readSimulationStatus(simulation)
     77     end
     78     filename = string(simulation.id, "/", simulation.id, ".$step.jld2")
     79     if verbose
     80         @info "Read simulation from $filename"
     81     end
     82     simulation = createSimulation()
     83     JLD2.@load(filename, simulation)
     84     return simulation
     85 end
     86 
     87 export writeSimulationStatus
     88 """
     89     writeSimulationStatus(simulation::Simulation;
     90                           folder::String=".",
     91                           verbose::Bool=false)
     92 
     93 Write current simulation status to disk in a minimal txt file.
     94 """
     95 function writeSimulationStatus(simulation::Simulation;
     96                                folder::String=".",
     97                                verbose::Bool=false)
     98     folder = folder * "/" * simulation.id
     99     mkpath(folder)
    100     filename = string(folder, "/", simulation.id, ".status.txt")
    101 
    102     writedlm(filename, [simulation.time
    103                         simulation.time/simulation.time_total*100.
    104                         float(simulation.file_number)])
    105     if verbose
    106         @info "Wrote status to $filename"
    107     end
    108     nothing
    109 end
    110 
    111 export readSimulationStatus
    112 """
    113     readSimulationStatus(simulation_id[, folder, verbose])
    114 
    115 Read the current simulation status from disk (`<sim.id>/<sim.id>.status.txt`)
    116 and return the last output file number.
    117 
    118 # Arguments
    119 * `simulation_id::String`: the simulation identifying string.
    120 * `folder::String="."`: the folder in which to search for the status file.
    121 * `verbose::Bool=true`: show simulation status in console.
    122 """
    123 function readSimulationStatus(simulation_id::String;
    124                               folder::String=".",
    125                               verbose::Bool=true)
    126 
    127     folder = folder * "/" * simulation_id
    128     filename = string(folder, "/", simulation_id, ".status.txt")
    129 
    130     data = readdlm(filename)
    131     if verbose
    132         @info "$simulation_id:\n" *
    133              "  time:             $(data[1]) s\n" *
    134              "  complete:         $(abs(data[2]))%\n" *
    135              "  last output file: $(Int(round(data[3])))\n"
    136     end
    137     return Int(round(data[3]))
    138 """
    139     readSimulationStatus(simulation[, folder, verbose])
    140 
    141 Read the current simulation status from disk (`<sim.id>/<sim.id>.status.txt`)
    142 and return the last output file number.
    143 
    144 # Arguments
    145 * `simulation::Simulation`: the simulation to read the status for.
    146 * `folder::String="."`: the folder in which to search for the status file.
    147 * `verbose::Bool=true`: show simulation status in console.
    148 """
    149 end
    150 function readSimulationStatus(sim::Simulation;
    151                               folder::String=".",
    152                               verbose::Bool=true)
    153     readSimulationStatus(sim.id, folder=folder, verbose=verbose)
    154 end
    155 
    156 export status
    157 """
    158     status(folder[, loop, t_int, colored_output, write_header, render)
    159 
    160 Shows the status of all simulations with output files written under the 
    161 specified `folder`, which is the current working directory by default.
    162 
    163 # Arguments
    164 `folder::String="."`: directory (including subdirectories) to scan for
    165     simulation output.
    166 `loop::Bool=false`: continue printing the status every `t_int` seconds.
    167 `t_int::Int=10`: interval between status updates when `loop=true`.
    168 `colored_output::Bool=true`: display output with colors.
    169 `write_header::Bool=true`: write header line explaining the data.
    170 visualize::Bool=false`: render the simulation output. Does not work well when
    171     `loop=true`, as the script regenerates (and overwrites)  all output graphics
    172     on every call.
    173 """
    174 function status(folder::String=".";
    175                 loop::Bool=false,
    176                 t_int::Int=10,
    177                 colored_output::Bool=true,
    178                 write_header::Bool=true,
    179                 visualize::Bool=false)
    180 
    181     if colored_output
    182         id_color_complete = :green
    183         id_color_in_progress = :yellow
    184         time_color = :white
    185         percentage_color = :blue
    186         lastfile_color = :cyan
    187     else
    188         id_color_complete = :default
    189         id_color_in_progress = :default
    190         time_color = :default
    191         percentage_color = :default
    192         lastfile_color = :default
    193     end
    194 
    195     repeat = true
    196     while repeat
    197 
    198         status_files = String[]
    199         println()
    200         println(Dates.format(Dates.DateTime(Dates.now()), Dates.RFC1123Format))
    201 
    202         for (root, dirs, files) in walkdir(folder, follow_symlinks=false)
    203 
    204             for file in files
    205                 if occursin(".status.txt", file)
    206                     push!(status_files, joinpath(root, file))
    207                 end
    208             end
    209         end
    210 
    211         if length(status_files) > 0
    212 
    213             cols = displaysize(stdout)[2]
    214             if write_header
    215                 for i=1:cols
    216                     print('-')
    217                 end
    218                 println()
    219 
    220                 left_label_width = 36
    221                 printstyled("simulation folder ", color=:default)
    222                 right_labels_width = 22
    223                 for i=18:cols-right_labels_width
    224                     print(' ')
    225                 end
    226                 printstyled("time  ", color=time_color)
    227                 printstyled("complete  ", color=percentage_color)
    228                 printstyled("files\n", color=lastfile_color)
    229                 for i=1:cols
    230                     print('-')
    231                 end
    232                 println()
    233             end
    234 
    235             for file in status_files
    236                 data = readdlm(file)
    237                 id = replace(file, ".status.txt" => "")
    238                 id = replace(id, "./" => "")
    239                 id = replace(id, r".*/" => "")
    240                 percentage = @sprintf "%9.0f%%" data[2]
    241                 lastfile = @sprintf "%7d" data[3]
    242                 if data[2] < 99.
    243                     printstyled("$id", color=id_color_in_progress)
    244                 else
    245                     printstyled("$id", color=id_color_complete)
    246                 end
    247                 right_fields_width = 25
    248                 for i=length(id):cols-right_fields_width
    249                     print(' ')
    250                 end
    251                 if data[1] < 60.0  # secs
    252                     time = @sprintf "%6.2fs" data[1]
    253                 elseif data[1] < 60.0*60.0  # mins
    254                     time = @sprintf "%6.2fm" data[1]/60.
    255                 elseif data[1] < 60.0*60.0*24.0  # hours
    256                     time = @sprintf "%6.2fh" data[1]/(60.0 * 60.0)
    257                 else  # days
    258                     time = @sprintf "%6.2fd" data[1]/(60.0 * 60.0 * 24.0)
    259                 end
    260                 printstyled("$time", color=time_color)
    261                 printstyled("$percentage", color=percentage_color)
    262                 printstyled("$lastfile\n", color=lastfile_color)
    263 
    264                 if visualize
    265                     sim = createSimulation(id)
    266                     render(sim)
    267                 end
    268             end
    269             if write_header
    270                 for i=1:cols
    271                     print('-')
    272                 end
    273                 println()
    274             end
    275         else
    276             @warn "no simulations found in $(pwd())/$folder"
    277         end
    278 
    279         if loop && t_int > 0
    280             sleep(t_int)
    281         end
    282         if !loop
    283             repeat = false
    284         end
    285     end
    286     nothing
    287 end
    288 
    289 export writeVTK
    290 """
    291 Write a VTK file to disk containing all grains in the `simulation` in an 
    292 unstructured mesh (file type `.vtu`).  These files can be read by ParaView and 
    293 can be visualized by applying a *Glyph* filter.
    294 
    295 If the simulation contains an `Ocean` data structure, it's contents will be 
    296 written to separate `.vtu` files.  This can be disabled by setting the argument 
    297 `ocean=false`.  The same is true for the atmosphere.
    298 
    299 The VTK files will be saved in a subfolder named after the simulation.
    300 """
    301 function writeVTK(simulation::Simulation;
    302                   folder::String=".",
    303                   verbose::Bool=true,
    304                   ocean::Bool=true,
    305                   atmosphere::Bool=true)
    306 
    307     simulation.file_number += 1
    308     folder = folder * "/" * simulation.id
    309     mkpath(folder)
    310 
    311     filename = string(folder, "/", simulation.id, ".grains.", 
    312                       simulation.file_number)
    313     writeGrainVTK(simulation, filename, verbose=verbose)
    314 
    315     filename = string(folder, "/", simulation.id, ".grain-interaction.", 
    316                       simulation.file_number)
    317     writeGrainInteractionVTK(simulation, filename, verbose=verbose)
    318 
    319     if typeof(simulation.ocean.input_file) != Bool && ocean
    320         filename = string(folder, "/", simulation.id, ".ocean.", 
    321                         simulation.file_number)
    322         writeGridVTK(simulation.ocean, filename, verbose=verbose)
    323     end
    324 
    325     if typeof(simulation.atmosphere.input_file) != Bool && atmosphere
    326         filename = string(folder, "/", simulation.id, ".atmosphere.", 
    327                         simulation.file_number)
    328         writeGridVTK(simulation.atmosphere, filename, verbose=verbose)
    329     end
    330     nothing
    331 end
    332 
    333 export writeGrainVTK
    334 """
    335 Write a VTK file to disk containing all grains in the `simulation` in an 
    336 unstructured mesh (file type `.vtu`).  These files can be read by ParaView and 
    337 can be visualized by applying a *Glyph* filter.  This function is called by 
    338 `writeVTK()`.
    339 """
    340 function writeGrainVTK(simulation::Simulation,
    341                          filename::String;
    342                          verbose::Bool=false)
    343 
    344     ifarr = convertGrainDataToArrays(simulation)
    345     
    346     # add arrays to VTK file
    347     vtkfile = WriteVTK.vtk_grid("$filename.vtu", ifarr.lin_pos,
    348                                 WriteVTK.MeshCell[])
    349 
    350     WriteVTK.vtk_point_data(vtkfile, ifarr.density, "Density [kg m^-3]")
    351 
    352     WriteVTK.vtk_point_data(vtkfile, ifarr.thickness, "Thickness [m]")
    353     WriteVTK.vtk_point_data(vtkfile, ifarr.contact_radius*2.,
    354                             "Diameter (contact) [m]")
    355     WriteVTK.vtk_point_data(vtkfile, ifarr.areal_radius*2.,
    356                             "Diameter (areal) [m]")
    357     WriteVTK.vtk_point_data(vtkfile, ifarr.circumreference,
    358                             "Circumreference  [m]")
    359     WriteVTK.vtk_point_data(vtkfile, ifarr.horizontal_surface_area,
    360                             "Horizontal surface area [m^2]")
    361     WriteVTK.vtk_point_data(vtkfile, ifarr.side_surface_area,
    362                             "Side surface area [m^2]")
    363     WriteVTK.vtk_point_data(vtkfile, ifarr.volume, "Volume [m^3]")
    364     WriteVTK.vtk_point_data(vtkfile, ifarr.mass, "Mass [kg]")
    365     WriteVTK.vtk_point_data(vtkfile, ifarr.moment_of_inertia,
    366                             "Moment of inertia [kg m^2]")
    367 
    368     WriteVTK.vtk_point_data(vtkfile, ifarr.lin_vel, "Linear velocity [m s^-1]")
    369     WriteVTK.vtk_point_data(vtkfile, ifarr.lin_acc,
    370                             "Linear acceleration [m s^-2]")
    371     WriteVTK.vtk_point_data(vtkfile, ifarr.force, "Sum of forces [N]")
    372     WriteVTK.vtk_point_data(vtkfile, ifarr.lin_disp, "Linear displacement [m]")
    373 
    374     WriteVTK.vtk_point_data(vtkfile, ifarr.ang_pos, "Angular position [rad]")
    375     WriteVTK.vtk_point_data(vtkfile, ifarr.ang_vel,
    376                             "Angular velocity [rad s^-1]")
    377     WriteVTK.vtk_point_data(vtkfile, ifarr.ang_acc,
    378                             "Angular acceleration [rad s^-2]")
    379     WriteVTK.vtk_point_data(vtkfile, ifarr.torque, "Sum of torques [N*m]")
    380 
    381     WriteVTK.vtk_point_data(vtkfile, ifarr.fixed, "Fixed in space [-]")
    382     WriteVTK.vtk_point_data(vtkfile, ifarr.allow_x_acc,
    383                             "Fixed but allow (x) acceleration [-]")
    384     WriteVTK.vtk_point_data(vtkfile, ifarr.allow_y_acc,
    385                             "Fixed but allow (y) acceleration [-]")
    386     WriteVTK.vtk_point_data(vtkfile, ifarr.allow_z_acc,
    387                             "Fixed but allow (z) acceleration [-]")
    388     WriteVTK.vtk_point_data(vtkfile, ifarr.rotating, "Free to rotate [-]")
    389     WriteVTK.vtk_point_data(vtkfile, ifarr.enabled, "Enabled [-]")
    390 
    391     WriteVTK.vtk_point_data(vtkfile, ifarr.contact_stiffness_normal,
    392                             "Contact stiffness (normal) [N m^-1]")
    393     WriteVTK.vtk_point_data(vtkfile, ifarr.contact_stiffness_tangential,
    394                             "Contact stiffness (tangential) [N m^-1]")
    395     WriteVTK.vtk_point_data(vtkfile, ifarr.contact_viscosity_normal,
    396                             "Contact viscosity (normal) [N m^-1 s]")
    397     WriteVTK.vtk_point_data(vtkfile, ifarr.contact_viscosity_tangential,
    398                             "Contact viscosity (tangential) [N m^-1 s]")
    399     WriteVTK.vtk_point_data(vtkfile, ifarr.contact_static_friction,
    400                             "Contact friction (static) [-]")
    401     WriteVTK.vtk_point_data(vtkfile, ifarr.contact_dynamic_friction,
    402                             "Contact friction (dynamic) [-]")
    403 
    404     WriteVTK.vtk_point_data(vtkfile, ifarr.youngs_modulus,
    405                             "Young's modulus [Pa]")
    406     WriteVTK.vtk_point_data(vtkfile, ifarr.poissons_ratio,
    407                             "Poisson's ratio [-]")
    408     WriteVTK.vtk_point_data(vtkfile, ifarr.tensile_strength,
    409                             "Tensile strength [Pa]")
    410     WriteVTK.vtk_point_data(vtkfile, ifarr.shear_strength,
    411                             "Shear strength [Pa]")
    412     WriteVTK.vtk_point_data(vtkfile, ifarr.strength_heal_rate,
    413                             "Bond strength healing rate [Pa/s]")
    414     WriteVTK.vtk_point_data(vtkfile, ifarr.fracture_toughness,
    415                             "Fracture toughness [m^0.5 Pa]")
    416 
    417     WriteVTK.vtk_point_data(vtkfile, ifarr.ocean_drag_coeff_vert,
    418                             "Ocean drag coefficient (vertical) [-]")
    419     WriteVTK.vtk_point_data(vtkfile, ifarr.ocean_drag_coeff_horiz,
    420                             "Ocean drag coefficient (horizontal) [-]")
    421     WriteVTK.vtk_point_data(vtkfile, ifarr.atmosphere_drag_coeff_vert,
    422                             "Atmosphere drag coefficient (vertical) [-]")
    423     WriteVTK.vtk_point_data(vtkfile, ifarr.atmosphere_drag_coeff_horiz,
    424                             "Atmosphere drag coefficient (horizontal) [-]")
    425 
    426     WriteVTK.vtk_point_data(vtkfile, ifarr.pressure,
    427                             "Contact pressure [Pa]")
    428     WriteVTK.vtk_point_data(vtkfile, ifarr.n_contacts,
    429                             "Number of contacts [-]")
    430 
    431     WriteVTK.vtk_point_data(vtkfile, ifarr.granular_stress,
    432                             "Granular stress [Pa]")
    433     WriteVTK.vtk_point_data(vtkfile, ifarr.ocean_stress,
    434                             "Ocean stress [Pa]")
    435     WriteVTK.vtk_point_data(vtkfile, ifarr.atmosphere_stress,
    436                             "Atmosphere stress [Pa]")
    437 
    438     WriteVTK.vtk_point_data(vtkfile, ifarr.thermal_energy,
    439                             "Thermal energy [J]")
    440 
    441     WriteVTK.vtk_point_data(vtkfile, ifarr.color,
    442                             "Color [-]")
    443 
    444     deleteGrainArrays!(ifarr)
    445     ifarr = 0
    446 
    447     outfiles = WriteVTK.vtk_save(vtkfile)
    448     if verbose
    449         @info "Output file: $(outfiles[1])"
    450     end
    451     nothing
    452 end
    453 
    454 export writeGrainInteractionVTK
    455 """
    456     writeGrainInteractionVTK(simulation::Simulation,
    457                                filename::String;
    458                                verbose::Bool=false)
    459 
    460 Saves grain interactions to `.vtp` files for visualization with VTK, for 
    461 example in Paraview.  Convert Cell Data to Point Data and use with Tube filter.
    462 """
    463 function writeGrainInteractionVTK(simulation::Simulation,
    464                                   filename::String;
    465                                   verbose::Bool=false)
    466 
    467     i1 = Int64[]
    468     i2 = Int64[]
    469     inter_particle_vector = Vector{Float64}[]
    470     force = Float64[]
    471     effective_radius = Float64[]
    472     contact_area = Float64[]
    473     contact_stiffness = Float64[]
    474     tensile_stress = Float64[]
    475     shear_displacement = Vector{Float64}[]
    476     contact_rotation = Vector{Float64}[]
    477     contact_age = Float64[]
    478     contact_area = Float64[]
    479     compressive_failure = Int[]
    480 
    481     for i=1:length(simulation.grains)
    482         for ic=1:simulation.Nc_max
    483             if simulation.grains[i].contacts[ic] > 0
    484                 j = simulation.grains[i].contacts[ic]
    485 
    486                 if !simulation.grains[i].enabled ||
    487                     !simulation.grains[j].enabled
    488                     continue
    489                 end
    490 
    491                 p = simulation.grains[i].lin_pos -
    492                     simulation.grains[j].lin_pos
    493                 dist = norm(p)
    494 
    495                 r_i = simulation.grains[i].contact_radius
    496                 r_j = simulation.grains[j].contact_radius
    497 
    498                 # skip visualization of contacts across periodic BCs
    499                 if dist > 5.0*(r_i + r_j) &&
    500                     (simulation.ocean.bc_west == 2 ||
    501                      simulation.ocean.bc_east == 2 ||
    502                      simulation.ocean.bc_north == 2 ||
    503                      simulation.ocean.bc_south == 2 ||
    504                      simulation.atmosphere.bc_west == 2 ||
    505                      simulation.atmosphere.bc_east == 2 ||
    506                      simulation.atmosphere.bc_north == 2 ||
    507                      simulation.atmosphere.bc_south == 2)
    508                     continue
    509                 end
    510                 δ_n = dist - (r_i + r_j)
    511                 R_ij = harmonicMean(r_i, r_j)
    512                 A_ij = R_ij*min(simulation.grains[i].thickness, 
    513                                 simulation.grains[j].thickness)
    514 
    515                 if simulation.grains[i].youngs_modulus > 0. &&
    516                     simulation.grains[j].youngs_modulus > 0.
    517                     E_ij = harmonicMean(simulation.grains[i].
    518                                         youngs_modulus,
    519                                         simulation.grains[j].
    520                                         youngs_modulus)
    521                     k_n = E_ij*A_ij/R_ij
    522                 else
    523                     k_n = harmonicMean(simulation.grains[i].
    524                                        contact_stiffness_normal,
    525                                        simulation.grains[j].
    526                                        contact_stiffness_normal)
    527                 end
    528 
    529                 
    530                 push!(i1, i)
    531                 push!(i2, j)
    532                 push!(inter_particle_vector, p)
    533 
    534                 push!(force, k_n*δ_n)
    535                 push!(effective_radius, R_ij)
    536                 push!(contact_area, A_ij)
    537                 push!(contact_stiffness, k_n)
    538                 push!(tensile_stress, k_n*δ_n/A_ij)
    539 
    540                 push!(shear_displacement,
    541                       simulation.grains[i]. contact_parallel_displacement[ic])
    542                 push!(contact_rotation,
    543                       simulation.grains[i].contact_rotation[ic])
    544 
    545                 push!(contact_age, simulation.grains[i].contact_age[ic])
    546                 push!(contact_area, simulation.grains[i].contact_area[ic])
    547                 push!(compressive_failure,
    548                       Int(simulation.grains[i].compressive_failure[ic]))
    549             end
    550         end
    551     end
    552 
    553     # Insert a piece for each grain interaction using grain positions as 
    554     # coordinates and connect them with lines by referencing their indexes.
    555     open(filename * ".vtp", "w") do f
    556         write(f, "<?xml version=\"1.0\"?>\n")
    557         write(f, "<VTKFile type=\"PolyData\" version=\"0.1\" " *
    558               "byte_order=\"LittleEndian\">\n")
    559         write(f, "  <PolyData>\n")
    560         write(f, "    <Piece " *
    561               "NumberOfPoints=\"$(length(simulation.grains))\" " *
    562               "NumberOfVerts=\"0\" " *
    563               "NumberOfLines=\"$(length(i1))\" " *
    564               "NumberOfStrips=\"0\" " *
    565               "NumberOfPolys=\"0\">\n")
    566         write(f, "      <PointData>\n")
    567         write(f, "      </PointData>\n")
    568         write(f, "      <CellData>\n")
    569 
    570         # Write values associated to each line
    571         write(f, "        <DataArray type=\"Float32\" " *
    572               "Name=\"Inter-particle vector [m]\" " *
    573               "NumberOfComponents=\"3\" format=\"ascii\">\n")
    574         for i=1:length(i1)
    575             write(f, "$(inter_particle_vector[i][1]) ")
    576             write(f, "$(inter_particle_vector[i][2]) ")
    577             write(f, "$(inter_particle_vector[i][3]) ")
    578         end
    579         write(f, "\n")
    580         write(f, "        </DataArray>\n")
    581         write(f, "        <DataArray type=\"Float32\" " *
    582               "Name=\"Shear displacement [m]\" " *
    583               "NumberOfComponents=\"3\" format=\"ascii\">\n")
    584         for i=1:length(i1)
    585             @inbounds write(f, "$(shear_displacement[i][1]) ")
    586             @inbounds write(f, "$(shear_displacement[i][2]) ")
    587             @inbounds write(f, "$(shear_displacement[i][3]) ")
    588         end
    589         write(f, "\n")
    590         write(f, "        </DataArray>\n")
    591 
    592         write(f, "        <DataArray type=\"Float32\" Name=\"Force [N]\" " *
    593               "NumberOfComponents=\"1\" format=\"ascii\">\n")
    594         for i=1:length(i1)
    595             @inbounds write(f, "$(force[i]) ")
    596         end
    597         write(f, "\n")
    598         write(f, "        </DataArray>\n")
    599 
    600         write(f, "        <DataArray type=\"Float32\" " *
    601               "Name=\"Effective radius [m]\" " *
    602               "NumberOfComponents=\"1\" format=\"ascii\">\n")
    603         for i=1:length(i1)
    604             @inbounds write(f, "$(effective_radius[i]) ")
    605         end
    606         write(f, "\n")
    607         write(f, "        </DataArray>\n")
    608 
    609         write(f, "        <DataArray type=\"Float32\" " *
    610               "Name=\"Contact area [m^2]\" " *
    611               "NumberOfComponents=\"1\" format=\"ascii\">\n")
    612         for i=1:length(i1)
    613             @inbounds write(f, "$(contact_area[i]) ")
    614         end
    615         write(f, "\n")
    616         write(f, "        </DataArray>\n")
    617 
    618         write(f, "        <DataArray type=\"Float32\" " *
    619               "Name=\"Contact stiffness [N/m]\" " *
    620               "NumberOfComponents=\"1\" format=\"ascii\">\n")
    621         for i=1:length(i1)
    622             @inbounds write(f, "$(contact_stiffness[i]) ")
    623         end
    624         write(f, "\n")
    625         write(f, "        </DataArray>\n")
    626 
    627         write(f, "        <DataArray type=\"Float32\" " *
    628               "Name=\"Tensile stress [Pa]\" " *
    629               "NumberOfComponents=\"1\" format=\"ascii\">\n")
    630         for i=1:length(i1)
    631             @inbounds write(f, "$(tensile_stress[i]) ")
    632         end
    633         write(f, "\n")
    634         write(f, "        </DataArray>\n")
    635 
    636         write(f, "        <DataArray type=\"Float32\" " *
    637               "Name=\"Contact rotation [rad]\" NumberOfComponents=\"3\" 
    638         format=\"ascii\">\n")
    639         for i=1:length(i1)
    640             @inbounds write(f, "$(contact_rotation[i][1]) ")
    641             @inbounds write(f, "$(contact_rotation[i][2]) ")
    642             @inbounds write(f, "$(contact_rotation[i][3]) ")
    643         end
    644         write(f, "\n")
    645         write(f, "        </DataArray>\n")
    646 
    647         write(f, "        <DataArray type=\"Float32\" " *
    648               "Name=\"Contact age [s]\" NumberOfComponents=\"1\" 
    649         format=\"ascii\">\n")
    650         for i=1:length(i1)
    651             @inbounds write(f, "$(contact_age[i]) ")
    652         end
    653         write(f, "\n")
    654         write(f, "        </DataArray>\n")
    655 
    656         write(f, "        <DataArray type=\"Float32\" " *
    657               "Name=\"Contact area [m*m]\" NumberOfComponents=\"1\" 
    658         format=\"ascii\">\n")
    659         for i=1:length(i1)
    660             @inbounds write(f, "$(contact_area[i]) ")
    661         end
    662         write(f, "\n")
    663         write(f, "        </DataArray>\n")
    664 
    665         write(f, "        <DataArray type=\"Float32\" " *
    666               "Name=\"Compressive failure [-]\" NumberOfComponents=\"1\" 
    667         format=\"ascii\">\n")
    668         for i=1:length(i1)
    669             @inbounds write(f, "$(Int(compressive_failure[i])) ")
    670         end
    671         write(f, "\n")
    672         write(f, "        </DataArray>\n")
    673 
    674         write(f, "      </CellData>\n")
    675         write(f, "      <Points>\n")
    676 
    677         # Write line endpoints (grain centers)
    678         #write(f, "        <DataArray Name=\"Position [m]\" type=\"Float32\" " *
    679         write(f, "        <DataArray type=\"Float32\" Name=\"Points\" " *
    680               "NumberOfComponents=\"3\" format=\"ascii\">\n")
    681         for i in simulation.grains
    682             @inbounds write(f, "$(i.lin_pos[1]) $(i.lin_pos[2]) $(i.lin_pos[3]) ")
    683         end
    684         write(f, "\n")
    685         write(f, "        </DataArray>\n")
    686         write(f, "      </Points>\n")
    687         write(f, "      <Verts>\n")
    688         write(f, "        <DataArray type=\"Int64\" Name=\"connectivity\" " *
    689               "format=\"ascii\">\n")
    690         write(f, "        </DataArray>\n")
    691         write(f, "        <DataArray type=\"Int64\" Name=\"offsets\" " *
    692               "format=\"ascii\">\n")
    693         write(f, "        </DataArray>\n")
    694         write(f, "      </Verts>\n")
    695         write(f, "      <Lines>\n")
    696 
    697         # Write contact connectivity by referring to point indexes
    698         write(f, "        <DataArray type=\"Int64\" Name=\"connectivity\" " *
    699               "format=\"ascii\">\n")
    700         for i=1:length(i1)
    701             @inbounds write(f, "$(i1[i] - 1) $(i2[i] - 1) ")
    702         end
    703         write(f, "\n")
    704         write(f, "        </DataArray>\n")
    705         
    706         # Write 0-indexed offset for the connectivity array for the end of each 
    707         # cell
    708         write(f, "        <DataArray type=\"Int64\" Name=\"offsets\" " *
    709               "format=\"ascii\">\n")
    710         for i=1:length(i1)
    711             @inbounds write(f, "$((i - 1)*2 + 2) ")
    712         end
    713         write(f, "\n")
    714         write(f, "        </DataArray>\n")
    715 
    716         write(f, "      </Lines>\n")
    717         write(f, "      <Strips>\n")
    718         write(f, "        <DataArray type=\"Int64\" Name=\"connectivity\" " *
    719               "format=\"ascii\">\n")
    720         write(f, "        </DataArray>\n")
    721         write(f, "        <DataArray type=\"Int64\" Name=\"offsets\" " *
    722               "format=\"ascii\">\n")
    723         write(f, "        </DataArray>\n")
    724         write(f, "      </Strips>\n")
    725         write(f, "      <Polys>\n")
    726         write(f, "        <DataArray type=\"Int64\" Name=\"connectivity\" " *
    727               "format=\"ascii\">\n")
    728         write(f, "        </DataArray>\n")
    729         write(f, "        <DataArray type=\"Int64\" Name=\"offsets\" " *
    730               "format=\"ascii\">\n")
    731         write(f, "        </DataArray>\n")
    732         write(f, "      </Polys>\n")
    733         write(f, "    </Piece>\n")
    734         write(f, "  </PolyData>\n")
    735         write(f, "</VTKFile>\n")
    736     end
    737     nothing
    738 end
    739 
    740 export writeOceanVTK
    741 """
    742 Write a VTK file to disk containing all ocean data in the `simulation` in a 
    743 structured grid (file type `.vts`).  These files can be read by ParaView and can 
    744 be visualized by applying a *Glyph* filter.  This function is called by 
    745 `writeVTK()`.
    746 """
    747 function writeGridVTK(grid::Any,
    748                       filename::String;
    749                       verbose::Bool=false)
    750     
    751     # make each coordinate array three-dimensional
    752     xq = similar(grid.u[:,:,:,1])
    753     yq = similar(grid.u[:,:,:,1])
    754     zq = similar(grid.u[:,:,:,1])
    755 
    756     for iz=1:size(xq, 3)
    757         @inbounds xq[:,:,iz] = grid.xq
    758         @inbounds yq[:,:,iz] = grid.yq
    759     end
    760     for ix=1:size(xq, 1)
    761         for iy=1:size(xq, 2)
    762             @inbounds zq[ix,iy,:] = grid.zl
    763         end
    764     end
    765 
    766     # add arrays to VTK file
    767     vtkfile = WriteVTK.vtk_grid("$filename.vts", xq, yq, zq)
    768 
    769     WriteVTK.vtk_point_data(vtkfile, grid.u[:, :, :, 1],
    770                             "u: Zonal velocity [m/s]")
    771     WriteVTK.vtk_point_data(vtkfile, grid.v[:, :, :, 1],
    772                             "v: Meridional velocity [m/s]")
    773     # write velocities as 3d vector
    774     vel = zeros(3, size(xq, 1), size(xq, 2), size(xq, 3))
    775     for ix=1:size(xq, 1)
    776         for iy=1:size(xq, 2)
    777             for iz=1:size(xq, 3)
    778                 @inbounds vel[1, ix, iy, iz] = grid.u[ix, iy, iz, 1]
    779                 @inbounds vel[2, ix, iy, iz] = grid.v[ix, iy, iz, 1]
    780             end
    781         end
    782     end
    783     WriteVTK.vtk_point_data(vtkfile, vel, "Velocity vector [m/s]")
    784 
    785     # Porosity is in the grids stored on the cell center, but is here
    786     # interpolated to the cell corners
    787     porosity = zeros(size(xq, 1), size(xq, 2), size(xq, 3))
    788     for ix=1:size(grid.xh, 1)
    789         for iy=1:size(grid.xh, 2)
    790             @inbounds porosity[ix, iy, 1] = grid.porosity[ix, iy]
    791             if ix == size(grid.xh, 1)
    792                 @inbounds porosity[ix+1, iy, 1] = grid.porosity[ix, iy]
    793             end
    794             if iy == size(grid.xh, 2)
    795                 @inbounds porosity[ix, iy+1, 1] = grid.porosity[ix, iy]
    796             end
    797             if ix == size(grid.xh, 1) && iy == size(grid.xh, 2)
    798                 @inbounds porosity[ix+1, iy+1, 1] = grid.porosity[ix, iy]
    799             end
    800         end
    801     end
    802     WriteVTK.vtk_point_data(vtkfile, porosity, "Porosity [-]")
    803 
    804     if typeof(grid) == Ocean
    805         WriteVTK.vtk_point_data(vtkfile, grid.h[:, :, :, 1],
    806                                 "h: Layer thickness [m]")
    807         WriteVTK.vtk_point_data(vtkfile, grid.e[:, :, :, 1],
    808                                 "e: Relative interface height [m]")
    809     end
    810 
    811     outfiles = WriteVTK.vtk_save(vtkfile)
    812     if verbose
    813         @info "Output file: $(outfiles[1])"
    814     end
    815     nothing
    816 end
    817 
    818 export writeParaviewPythonScript
    819 """
    820 function writeParaviewPythonScript(simulation,
    821                                    [filename, folder, vtk_folder, verbose])
    822 
    823 Create a `".py"` script for visualizing the simulation VTK files in Paraview.
    824 The script can be run from the command line with `pvpython` (bundled with
    825 Paraview), or from the interactive Python shell inside Paraview.
    826 
    827 # Arguments
    828 * `simulation::Simulation`: input simulation file containing the data.
    829 * `filename::String`: output file name for the Python script. At its default
    830     (blank) value, the script is named after the simulation id (`simulation.id`).
    831 * `folder::String`: output directory, current directory the default.
    832 * `vtk_folder::String`: directory containing the VTK output files, by default
    833     points to the full system path equivalent to `"./<simulation.id>/"`.
    834 * `save_animation::Bool`: make the generated script immediately save a rendered
    835     animation to disk when the `".py"` script is called.
    836 * `verbose::Bool`: show diagnostic information during
    837 function call, on by
    838     default.
    839 """
    840 function writeParaviewPythonScript(simulation::Simulation;
    841                                    filename::String="",
    842                                    folder::String=".",
    843                                    vtk_folder::String="",
    844                                    save_animation::Bool=true,
    845                                    save_images::Bool=false,
    846                                    width::Integer=1920,
    847                                    height::Integer=1080,
    848                                    framerate::Integer=10,
    849                                    grains_color_scheme::String="X Ray",
    850                                    verbose::Bool=true)
    851     if filename == ""
    852         folder = string(folder, "/", simulation.id)
    853         mkpath(folder)
    854         filename = string(folder, "/", simulation.id, ".py")
    855     end
    856     if vtk_folder == ""
    857         vtk_folder = string(pwd(), "/", simulation.id)
    858     end
    859 
    860     if simulation.file_number == 0
    861         simulation.file_number = readSimulationStatus(simulation,
    862                                                       verbose=verbose)
    863     end
    864 
    865     open(filename, "w") do f
    866         write(f, """from paraview.simple import *
    867 paraview.simple._DisableFirstRenderCameraReset()
    868 FileName=[""")
    869         for i=1:simulation.file_number
    870             write(f, "'$(vtk_folder)/$(simulation.id).grains.$(i).vtu', ")
    871         end
    872         write(f, """]
    873 imagegrains = XMLUnstructuredGridReader(FileName=FileName)
    874 
    875 imagegrains.PointArrayStatus = [
    876 'Density [kg m^-3]',
    877 'Thickness [m]',
    878 'Diameter (contact) [m]',
    879 'Diameter (areal) [m]',
    880 'Circumreference  [m]',
    881 'Horizontal surface area [m^2]',
    882 'Side surface area [m^2]',
    883 'Volume [m^3]',
    884 'Mass [kg]',
    885 'Moment of inertia [kg m^2]',
    886 'Linear velocity [m s^-1]',
    887 'Linear acceleration [m s^-2]',
    888 'Sum of forces [N]',
    889 'Linear displacement [m]',
    890 'Angular position [rad]',
    891 'Angular velocity [rad s^-1]',
    892 'Angular acceleration [rad s^-2]',
    893 'Sum of torques [N*m]',
    894 'Fixed in space [-]',
    895 'Fixed but allow (x) acceleration [-]',
    896 'Fixed but allow (y) acceleration [-]',
    897 'Fixed but allow (z) acceleration [-]',
    898 'Free to rotate [-]',
    899 'Enabled [-]',
    900 'Contact stiffness (normal) [N m^-1]',
    901 'Contact stiffness (tangential) [N m^-1]',
    902 'Contact viscosity (normal) [N m^-1 s]',
    903 'Contact viscosity (tangential) [N m^-1 s]',
    904 'Contact friction (static) [-]',
    905 'Contact friction (dynamic) [-]',
    906 "Young's modulus [Pa]",
    907 "Poisson's ratio [-]",
    908 'Tensile strength [Pa]'
    909 'Shear strength [Pa]'
    910 'Strength heal rate [Pa/s]'
    911 'Compressive strength prefactor [m^0.5 Pa]',
    912 'Ocean drag coefficient (vertical) [-]',
    913 'Ocean drag coefficient (horizontal) [-]',
    914 'Atmosphere drag coefficient (vertical) [-]',
    915 'Atmosphere drag coefficient (horizontal) [-]',
    916 'Contact pressure [Pa]',
    917 'Number of contacts [-]',
    918 'Granular stress [Pa]',
    919 'Ocean stress [Pa]',
    920 'Atmosphere stress [Pa]',
    921 'Color [-]']
    922 
    923 animationScene1 = GetAnimationScene()
    924 
    925 # update animation scene based on data timesteps
    926 animationScene1.UpdateAnimationUsingDataTimeSteps()
    927 
    928 # get active view
    929 renderView1 = GetActiveViewOrCreate('RenderView')
    930 # uncomment following to set a specific view size
    931 # renderView1.ViewSize = [2478, 1570]
    932 
    933 # show data in view
    934 imagegrainsDisplay = Show(imagegrains, renderView1)
    935 # trace defaults for the display properties.
    936 imagegrainsDisplay.Representation = 'Surface'
    937 imagegrainsDisplay.AmbientColor = [0.0, 0.0, 0.0]
    938 imagegrainsDisplay.ColorArrayName = [None, '']
    939 imagegrainsDisplay.OSPRayScaleArray = 'Angular acceleration [rad s^-2]'
    940 imagegrainsDisplay.OSPRayScaleFunction = 'PiecewiseFunction'
    941 imagegrainsDisplay.SelectOrientationVectors = 'Angular acceleration [rad s^-2]'
    942 imagegrainsDisplay.ScaleFactor = 6.050000000000001
    943 imagegrainsDisplay.SelectScaleArray = 'Angular acceleration [rad s^-2]'
    944 imagegrainsDisplay.GlyphType = 'Arrow'
    945 imagegrainsDisplay.GlyphTableIndexArray = 'Angular acceleration [rad s^-2]'
    946 imagegrainsDisplay.DataAxesGrid = 'GridAxesRepresentation'
    947 imagegrainsDisplay.PolarAxes = 'PolarAxesRepresentation'
    948 imagegrainsDisplay.ScalarOpacityUnitDistance = 64.20669746996803
    949 imagegrainsDisplay.GaussianRadius = 3.0250000000000004
    950 imagegrainsDisplay.SetScaleArray = ['POINTS', 'Atmosphere drag coefficient (horizontal) [-]']
    951 imagegrainsDisplay.ScaleTransferFunction = 'PiecewiseFunction'
    952 imagegrainsDisplay.OpacityArray = ['POINTS', 'Atmosphere drag coefficient (horizontal) [-]']
    953 imagegrainsDisplay.OpacityTransferFunction = 'PiecewiseFunction'
    954 
    955 # init the 'GridAxesRepresentation' selected for 'DataAxesGrid'
    956 imagegrainsDisplay.DataAxesGrid.XTitleColor = [0.0, 0.0, 0.0]
    957 imagegrainsDisplay.DataAxesGrid.YTitleColor = [0.0, 0.0, 0.0]
    958 imagegrainsDisplay.DataAxesGrid.ZTitleColor = [0.0, 0.0, 0.0]
    959 imagegrainsDisplay.DataAxesGrid.GridColor = [0.0, 0.0, 0.0]
    960 imagegrainsDisplay.DataAxesGrid.XLabelColor = [0.0, 0.0, 0.0]
    961 imagegrainsDisplay.DataAxesGrid.YLabelColor = [0.0, 0.0, 0.0]
    962 imagegrainsDisplay.DataAxesGrid.ZLabelColor = [0.0, 0.0, 0.0]
    963 
    964 # init the 'PolarAxesRepresentation' selected for 'PolarAxes'
    965 imagegrainsDisplay.PolarAxes.PolarAxisTitleColor = [0.0, 0.0, 0.0]
    966 imagegrainsDisplay.PolarAxes.PolarAxisLabelColor = [0.0, 0.0, 0.0]
    967 imagegrainsDisplay.PolarAxes.LastRadialAxisTextColor = [0.0, 0.0, 0.0]
    968 imagegrainsDisplay.PolarAxes.SecondaryRadialAxesTextColor = [0.0, 0.0, 0.0]
    969 
    970 # reset view to fit data
    971 renderView1.ResetCamera()
    972 
    973 #changing interaction mode based on data extents
    974 renderView1.InteractionMode = '2D'
    975 
    976 # update the view to ensure updated data information
    977 renderView1.Update()
    978 
    979 # create a new 'Glyph'
    980 glyph1 = Glyph(Input=imagegrains,
    981     GlyphType='Arrow')
    982 glyph1.add_attribute = ['POINTS', 'Atmosphere drag coefficient (horizontal) [-]']
    983 glyph1.add_attribute = ['POINTS', 'Angular acceleration [rad s^-2]']
    984 glyph1.ScaleFactor = 6.050000000000001
    985 glyph1.GlyphTransform = 'Transform2'
    986 
    987 # Properties modified on glyph1
    988 glyph1.add_attribute = ['POINTS', 'Diameter (areal) [m]']
    989 glyph1.add_attribute = ['POINTS', 'Angular position [rad]']
    990 glyph1.ScaleFactor = 1.0
    991 glyph1.GlyphMode = 'All Points'
    992 
    993 # get color transfer function/color map for 'Diameterarealm'
    994 diameterarealmLUT = GetColorTransferFunction('Diameterarealm')
    995 
    996 # show data in view
    997 glyph1Display = Show(glyph1, renderView1)
    998 # trace defaults for the display properties.
    999 glyph1Display.Representation = 'Surface'
   1000 glyph1Display.AmbientColor = [0.0, 0.0, 0.0]
   1001 glyph1Display.ColorArrayName = ['POINTS', 'Diameter (areal) [m]']
   1002 glyph1Display.LookupTable = diameterarealmLUT
   1003 glyph1Display.OSPRayScaleArray = 'Diameter (areal) [m]'
   1004 glyph1Display.OSPRayScaleFunction = 'PiecewiseFunction'
   1005 glyph1Display.SelectOrientationVectors = 'GlyphVector'
   1006 glyph1Display.ScaleFactor = 6.1000000000000005
   1007 glyph1Display.SelectScaleArray = 'Diameter (areal) [m]'
   1008 glyph1Display.GlyphType = 'Arrow'
   1009 glyph1Display.GlyphTableIndexArray = 'Diameter (areal) [m]'
   1010 glyph1Display.DataAxesGrid = 'GridAxesRepresentation'
   1011 glyph1Display.PolarAxes = 'PolarAxesRepresentation'
   1012 glyph1Display.GaussianRadius = 3.0500000000000003
   1013 glyph1Display.SetScaleArray = ['POINTS', 'Diameter (areal) [m]']
   1014 glyph1Display.ScaleTransferFunction = 'PiecewiseFunction'
   1015 glyph1Display.OpacityArray = ['POINTS', 'Diameter (areal) [m]']
   1016 glyph1Display.OpacityTransferFunction = 'PiecewiseFunction'
   1017 
   1018 # init the 'GridAxesRepresentation' selected for 'DataAxesGrid'
   1019 glyph1Display.DataAxesGrid.XTitleColor = [0.0, 0.0, 0.0]
   1020 glyph1Display.DataAxesGrid.YTitleColor = [0.0, 0.0, 0.0]
   1021 glyph1Display.DataAxesGrid.ZTitleColor = [0.0, 0.0, 0.0]
   1022 glyph1Display.DataAxesGrid.GridColor = [0.0, 0.0, 0.0]
   1023 glyph1Display.DataAxesGrid.XLabelColor = [0.0, 0.0, 0.0]
   1024 glyph1Display.DataAxesGrid.YLabelColor = [0.0, 0.0, 0.0]
   1025 glyph1Display.DataAxesGrid.ZLabelColor = [0.0, 0.0, 0.0]
   1026 
   1027 # init the 'PolarAxesRepresentation' selected for 'PolarAxes'
   1028 glyph1Display.PolarAxes.PolarAxisTitleColor = [0.0, 0.0, 0.0]
   1029 glyph1Display.PolarAxes.PolarAxisLabelColor = [0.0, 0.0, 0.0]
   1030 glyph1Display.PolarAxes.LastRadialAxisTextColor = [0.0, 0.0, 0.0]
   1031 glyph1Display.PolarAxes.SecondaryRadialAxesTextColor = [0.0, 0.0, 0.0]
   1032 
   1033 # show color bar/color legend
   1034 glyph1Display.SetScalarBarVisibility(renderView1, True)
   1035 
   1036 # update the view to ensure updated data information
   1037 renderView1.Update()
   1038 
   1039 # reset view to fit data
   1040 renderView1.ResetCamera()
   1041 
   1042 # Properties modified on glyph1
   1043 glyph1.GlyphType = 'Sphere'
   1044 
   1045 # update the view to ensure updated data information
   1046 renderView1.Update()
   1047 
   1048 # hide color bar/color legend
   1049 glyph1Display.SetScalarBarVisibility(renderView1, False)
   1050 
   1051 # rescale color and/or opacity maps used to exactly fit the current data range
   1052 glyph1Display.RescaleTransferFunctionToDataRange(False, True)
   1053 
   1054 # get opacity transfer function/opacity map for 'Diameterarealm'
   1055 diameterarealmPWF = GetOpacityTransferFunction('Diameterarealm')
   1056 
   1057 # Apply a preset using its name. Note this may not work as expected when presets have duplicate names.
   1058 diameterarealmLUT.ApplyPreset('$(grains_color_scheme)', True)
   1059 
   1060 # Hide orientation axes
   1061 renderView1.OrientationAxesVisibility = 0
   1062 
   1063 # current camera placement for renderView1
   1064 renderView1.InteractionMode = '2D'
   1065 """)
   1066         if save_animation
   1067             write(f, """
   1068 SaveAnimation('$(vtk_folder)/$(simulation.id).avi', renderView1,
   1069 ImageResolution=[$(width), $(height)],
   1070 FrameRate=$(framerate),
   1071 FrameWindow=[0, $(simulation.file_number)])
   1072 """)
   1073         end
   1074 
   1075         if save_images
   1076             write(f, """
   1077 SaveAnimation('$(folder)/$(simulation.id).png', renderView1,
   1078 ImageResolution=[$(width), $(height)],
   1079 FrameRate=$(framerate),
   1080 FrameWindow=[0, $(simulation.file_number)])
   1081 """)
   1082         end
   1083     end
   1084     if verbose
   1085         @info "$(filename) written, execute with " *
   1086              "'pvpython $(vtk_folder)/$(simulation.id).py'"
   1087     end
   1088 end
   1089 
   1090 export render
   1091 """
   1092     render(simulation[, pvpython, images, animation])
   1093 
   1094 Wrapper function which calls `writeParaviewPythonScript(...)` and executes it
   1095 from the shell using the supplied `pvpython` argument.
   1096 
   1097 # Arguments
   1098 * `simulation::Simulation`: simulation object containing the grain data.
   1099 * `pvpython::String`: path to the `pvpython` executable to use.  By default, the
   1100     script uses the pvpython in the system PATH.
   1101 * `images::Bool`: render images to disk (default: true)
   1102 * `gif::Bool`: merge images as GIF and save to disk (default: false, requires
   1103     `images=true`)
   1104 * `animation::Bool`: render animation as movie to disk (default: false). If
   1105     ffmpeg is available on the system, the `.avi` file is converted to `.mp4`.
   1106 * `trim::Bool`: trim images in animated sequence (default: true)
   1107 * `reverse::Bool`: if `images=true` additionally render reverse-animated gif
   1108     (default: false)
   1109 """
   1110 function render(simulation::Simulation; pvpython::String="pvpython",
   1111                 images::Bool=true,
   1112                 gif::Bool=false,
   1113                 animation::Bool=false,
   1114                 trim::Bool=true,
   1115                 reverse::Bool=false)
   1116 
   1117     writeParaviewPythonScript(simulation, save_animation=animation,
   1118                               save_images=images, verbose=false)
   1119     try
   1120         run(`$(pvpython) $(simulation.id)/$(simulation.id).py`)
   1121 
   1122         if animation
   1123             try
   1124                 run(`ffmpeg -i $(simulation.id)/$(simulation.id).avi 
   1125                     -vf scale='trunc\(iw/2\)\*2:trunc\(ih/2\)\*2' 
   1126                     -c:v libx264 -profile:v high -pix_fmt yuv420p 
   1127                     -g 30 -r 30 -y 
   1128                     $(simulation.id)/$(simulation.id).mp4`)
   1129                 if isfile("$(simulation.id)/$(simulation.id).mp4")
   1130                     rm("$(simulation.id)/$(simulation.id).avi")
   1131                 end
   1132             catch return_signal
   1133                 if isa(return_signal, Base.IOError)
   1134                     @warn "Could not run external ffmpeg command, " *
   1135                     "skipping conversion from " *
   1136                     "$(simulation.id)/$(simulation.id).avi to mp4."
   1137                 end
   1138             end
   1139         end
   1140 
   1141         # if available, use imagemagick to create gif from images
   1142         if images && gif
   1143             try
   1144                 trim_string = ""
   1145                 if trim
   1146                     trim_string = "-trim"
   1147                 end
   1148 
   1149                 # use ImageMagick installed with Homebrew.jl if available,
   1150                 # otherwise search for convert in $PATH
   1151                 convert = "convert"
   1152 
   1153                 run(`$convert $trim_string +repage -delay 10 
   1154                     -transparent-color white 
   1155                     -loop 0 $(simulation.id)/$(simulation.id).'*'.png 
   1156                     $(simulation.id)/$(simulation.id).gif`)
   1157                 if reverse
   1158                     run(`$convert -trim +repage -delay 10 
   1159                         -transparent-color white 
   1160                         -loop 0 -reverse
   1161                         $(simulation.id)/$(simulation.id).'*'.png 
   1162                         $(simulation.id)/$(simulation.id)-reverse.gif`)
   1163                 end
   1164             catch return_signal
   1165                 if isa(return_signal, Base.IOError)
   1166                     @warn "Skipping gif merge since `$convert` " *
   1167                         "was not found."
   1168                 end
   1169             end
   1170         end
   1171     catch return_signal
   1172         if isa(return_signal, Base.IOError)
   1173             error("`pvpython` was not found.")
   1174         end
   1175     end
   1176 end
   1177 
   1178 export removeSimulationFiles
   1179 """
   1180     removeSimulationFiles(simulation[, folder])
   1181 
   1182 Remove all simulation output files from the specified folder.
   1183 """
   1184 function removeSimulationFiles(simulation::Simulation; folder::String=".")
   1185     folder = folder * "/" * simulation.id
   1186     run(`sh -c "rm -rf $(folder)/$(simulation.id).*.vtu"`)
   1187     run(`sh -c "rm -rf $(folder)/$(simulation.id).*.vtp"`)
   1188     run(`sh -c "rm -rf $(folder)/$(simulation.id).*.vts"`)
   1189     run(`sh -c "rm -rf $(folder)/$(simulation.id).status.txt"`)
   1190     run(`sh -c "rm -rf $(folder)/$(simulation.id).*.jld2"`)
   1191     run(`sh -c "rm -rf $(folder)/$(simulation.id).py"`)
   1192     run(`sh -c "rm -rf $(folder)/$(simulation.id).avi"`)
   1193     run(`sh -c "rm -rf $(folder)/$(simulation.id).*.png"`)
   1194     run(`sh -c "rm -rf $(folder)"`)
   1195     nothing
   1196 end
   1197 
   1198 export plotGrainSizeDistribution
   1199 """
   1200     plotGrainSizeDistribution(simulation, [filename_postfix, nbins,
   1201                                 size_type, filetype, gnuplot_terminal,
   1202                                 skip_fixed, log_y, verbose)
   1203 
   1204 Plot the grain size distribution as a histogram and save it to the disk.  The 
   1205 plot is saved accoring to the simulation id, the optional `filename_postfix` 
   1206 string, and the `filetype`, and is written to the current folder.
   1207 
   1208 # Arguments
   1209 * `simulation::Simulation`: the simulation object containing the grains.
   1210 * `filename_postfix::String`: optional string for the output filename.
   1211 * `nbins::Int`: number of bins in the histogram (default = 12).
   1212 * `size_type::String`: specify whether to use the `contact` or `areal` radius 
   1213     for the grain size.  The default is `contact`.
   1214 * `filetype::String`: the output file type (default = "png").
   1215 * `gnuplot_terminal::String`: the gnuplot output terminal to use (default =
   1216     "png").
   1217 * `skip_fixed::Bool`: ommit grains that are fixed in space from the size 
   1218     distribution (default = true).
   1219 * `log_y::Bool`: plot y-axis in log scale.
   1220 * `verbose::String`: show output file as info message in stdout (default = 
   1221     true).
   1222 """
   1223 function plotGrainSizeDistribution(simulation::Simulation;
   1224                                      filename_postfix::String = "",
   1225                                      nbins::Int=12,
   1226                                      size_type::String = "contact",
   1227                                      filetype::String = "png",
   1228                                      gnuplot_terminal::String = "png",
   1229                                      skip_fixed::Bool = true,
   1230                                      log_y::Bool = false,
   1231                                      verbose::Bool = true)
   1232 
   1233     diameters = Float64[]
   1234     for i=1:length(simulation.grains)
   1235         if simulation.grains[i].fixed && skip_fixed
   1236             continue
   1237         end
   1238         if size_type == "contact"
   1239             push!(diameters, simulation.grains[i].contact_radius*2.)
   1240         elseif size_type == "areal"
   1241             push!(diameters, simulation.grains[i].areal_radius*2.)
   1242         else
   1243             error("size_type '$size_type' not understood")
   1244         end
   1245     end
   1246 
   1247     filename = string(simulation.id * filename_postfix * 
   1248                       "-grain-size-distribution." * filetype)
   1249 
   1250     # write data to temporary file on disk
   1251     datafile = Base.Filesystem.tempname()
   1252     writedlm(datafile, diameters)
   1253     gnuplotscript = Base.Filesystem.tempname()
   1254 
   1255     #if maximum(diameters) ≈ minimum(diameters)
   1256         #@info "Overriding `nbins = $nbins` -> `nbins = 1`."
   1257         #nbins = 1
   1258     #end
   1259 
   1260     open(gnuplotscript, "w") do f
   1261 
   1262         write(f, """#!/usr/bin/env gnuplot
   1263               set term $gnuplot_terminal
   1264               set out "$(filename)"\n""")
   1265         if log_y
   1266             write(f, "set logscale y\n")
   1267         end
   1268         write(f, """set xlabel "Diameter [m]"
   1269               set ylabel "Count [-]"
   1270               binwidth = $((maximum(diameters) - minimum(diameters)+1e-7)/nbins)
   1271               binstart = $(minimum(diameters))
   1272               set boxwidth 1.0*binwidth
   1273               set style fill solid 0.5
   1274               set key off
   1275               hist = 'u (binwidth*(floor((\$1-binstart)/binwidth)+0.5)+binstart):(1.0) smooth freq w boxes'
   1276               plot "$(datafile)" i 0 @hist ls 1
   1277               """)
   1278     end
   1279 
   1280     try
   1281         run(`gnuplot $gnuplotscript`)
   1282     catch return_signal
   1283         if isa(return_signal, Base.IOError)
   1284             error("Could not launch external gnuplot process")
   1285         end
   1286     end
   1287 
   1288     if verbose
   1289         @info filename
   1290     end
   1291 end
   1292 
   1293 export plotGrains
   1294 """
   1295     plotGrains(simulation, [filetype, gnuplot_terminal, verbose])
   1296 
   1297 Plot the grains using Gnuplot and save the figure to disk.
   1298 
   1299 # Arguments
   1300 * `simulation::Simulation`: the simulation object containing the grains.
   1301 * `filetype::String`: the output file type (default = "png").
   1302 * `gnuplot_terminal::String`: the gnuplot output terminal to use (default =
   1303     "png crop size 1200,1200").
   1304 * `plot_interactions::Bool`: show grain-grain interactions in the plot.
   1305 * `verbose::String`: show output file as info message in stdout (default = 
   1306     true).
   1307 """
   1308 function plotGrains(sim::Simulation;
   1309                     filetype::String = "png",
   1310                     gnuplot_terminal::String = "png crop size 1200,1200",
   1311                     plot_interactions::Bool = true,
   1312                     palette_scalar::String = "contact_radius",
   1313                     cbrange::Vector{Float64} = [NaN, NaN],
   1314                     show_figure::Bool = true,
   1315                     verbose::Bool = true)
   1316 
   1317     mkpath(sim.id)
   1318     filename = string(sim.id, "/", sim.id, ".grains.", sim.file_number, ".",
   1319                       filetype)
   1320 
   1321     # prepare grain data
   1322     x = Float64[]
   1323     y = Float64[]
   1324     r = Float64[]
   1325     scalars = Float64[]
   1326     for grain in sim.grains
   1327         push!(x, grain.lin_pos[1])
   1328         push!(y, grain.lin_pos[2])
   1329         push!(r, grain.contact_radius)
   1330 
   1331         if palette_scalar == "contact_radius"
   1332             push!(scalars, grain.contact_radius)
   1333 
   1334         elseif palette_scalar == "areal_radius"
   1335             push!(scalars, grain.areal_radius)
   1336 
   1337         elseif palette_scalar == "color"
   1338             push!(scalars, grain.color)
   1339 
   1340         else
   1341             error("palette_scalar = '$palette_scalar' not understood.")
   1342         end
   1343     end
   1344 
   1345     # prepare interaction data
   1346     if plot_interactions
   1347         i1 = Int64[]
   1348         i2 = Int64[]
   1349         inter_particle_vector = Vector{Float64}[]
   1350         force = Float64[]
   1351         effective_radius = Float64[]
   1352         contact_area = Float64[]
   1353         contact_stiffness = Float64[]
   1354         tensile_stress = Float64[]
   1355         shear_displacement = Vector{Float64}[]
   1356         contact_rotation = Vector{Float64}[]
   1357         contact_age = Float64[]
   1358         compressive_failure = Int[]
   1359         for i=1:length(sim.grains)
   1360             for ic=1:sim.Nc_max
   1361                 if sim.grains[i].contacts[ic] > 0
   1362                     j = sim.grains[i].contacts[ic]
   1363 
   1364                     if !sim.grains[i].enabled ||
   1365                         !sim.grains[j].enabled
   1366                         continue
   1367                     end
   1368 
   1369                     p = sim.grains[i].lin_pos -
   1370                         sim.grains[j].lin_pos
   1371                     dist = norm(p)
   1372 
   1373                     r_i = sim.grains[i].contact_radius
   1374                     r_j = sim.grains[j].contact_radius
   1375                     δ_n = dist - (r_i + r_j)
   1376                     R_ij = harmonicMean(r_i, r_j)
   1377 
   1378                     if sim.grains[i].youngs_modulus > 0. &&
   1379                         sim.grains[j].youngs_modulus > 0.
   1380                         E_ij = harmonicMean(sim.grains[i].
   1381                                             youngs_modulus,
   1382                                             sim.grains[j].
   1383                                             youngs_modulus)
   1384                         A_ij = R_ij*min(sim.grains[i].thickness, 
   1385                                         sim.grains[j].thickness)
   1386                         k_n = E_ij*A_ij/R_ij
   1387                     else
   1388                         k_n = harmonicMean(sim.grains[i].
   1389                                            contact_stiffness_normal,
   1390                                            sim.grains[j].
   1391                                            contact_stiffness_normal)
   1392                     end
   1393 
   1394                     push!(i1, i)
   1395                     push!(i2, j)
   1396                     push!(inter_particle_vector, p)
   1397 
   1398                     push!(force, k_n*δ_n)
   1399                     push!(effective_radius, R_ij)
   1400                     push!(contact_area, A_ij)
   1401                     push!(contact_stiffness, k_n)
   1402                     push!(tensile_stress, k_n*δ_n/A_ij)
   1403 
   1404                     push!(shear_displacement,
   1405                           sim.grains[i].contact_parallel_displacement[ic])
   1406                     push!(contact_rotation,
   1407                           sim.grains[i].contact_rotation[ic])
   1408 
   1409                     push!(contact_age, sim.grains[i].contact_age[ic])
   1410                     push!(compressive_failure,
   1411                       sim.grains[i].compressive_failure[ic])
   1412                 end
   1413             end
   1414         end
   1415     end
   1416 
   1417     # write grain data to temporary file on disk
   1418     datafile = Base.Filesystem.tempname()
   1419     writedlm(datafile, [x y r scalars])
   1420     gnuplotscript = Base.Filesystem.tempname()
   1421 
   1422     #=
   1423     if plot_interactions
   1424         # write grain data to temporary file on disk
   1425         datafile_int = Base.Filesystem.tempname()
   1426 
   1427         open(datafile_int, "w") do f
   1428             for i=1:length(i1)
   1429                 write(f, "$(sim.grains[i1[i]].lin_pos[1]) ")
   1430                 write(f, "$(sim.grains[i1[i]].lin_pos[2]) ")
   1431                 write(f, "$(sim.grains[i2[i]].lin_pos[1]) ")
   1432                 write(f, "$(sim.grains[i2[i]].lin_pos[2]) ")
   1433                 write(f, "$(force[i]) ")
   1434                 write(f, "$(effective_radius[i]) ")
   1435                 write(f, "$(contact_area[i]) ")
   1436                 write(f, "$(contact_stiffness[i]) ")
   1437                 write(f, "$(tensile_stress[i]) ")
   1438                 write(f, "$(shear_displacement[i]) ")
   1439                 write(f, "$(contact_age[i]) ")
   1440                 write(f, "\n")
   1441             end
   1442         end
   1443     end=#
   1444 
   1445     open(gnuplotscript, "w") do f
   1446 
   1447         write(f, """#!/usr/bin/env gnuplot
   1448               set term $(gnuplot_terminal)
   1449               set out "$(filename)"
   1450               set xlabel "x [m]"
   1451               set ylabel "y [m]"\n""")
   1452         if typeof(sim.ocean.input_file) != Bool
   1453             write(f, "set xrange ")
   1454             write(f, "[$(sim.ocean.xq[1,1]):$(sim.ocean.xq[end,end])]\n")
   1455             write(f, "set yrange ")
   1456             write(f, "[$(sim.ocean.yq[1,1]):$(sim.ocean.yq[end,end])]\n")
   1457         else
   1458             write(f, "set xrange [$(minimum(x - r)):$(maximum(x + r))]\n")
   1459             write(f, "set yrange [$(minimum(y - r)):$(maximum(y + r))]\n")
   1460         end
   1461 
   1462         # light gray to black to red
   1463         #write(f, "set palette defined ( 1 '#d3d3d3', 2 '#000000')\n")
   1464 
   1465         # light gray to black to red
   1466         write(f, "set palette defined ( 1 '#d3d3d3', 2 '#000000', 3 '#993333')\n")
   1467         
   1468         if !isnan(cbrange[1])
   1469             write(f, "set cbrange [$(cbrange[1]):$(cbrange[2])]\n")
   1470         end
   1471 
   1472         # gray to white
   1473         #write(f, "set palette defined (0 'gray', 1 'white')\n")
   1474 
   1475         write(f, """set cblabel "$palette_scalar"
   1476               set size ratio -1
   1477               set key off\n""")
   1478 
   1479         if length(i1) > 0
   1480             max_tensile_stress = maximum(abs.(tensile_stress))
   1481             max_line_width = 5.
   1482             if plot_interactions
   1483                 write(f, "set cbrange [-$max_tensile_stress:$max_tensile_stress]\n")
   1484                 for i=1:length(i1)
   1485                     write(f, "set arrow $i from " *
   1486                           "$(sim.grains[i1[i]].lin_pos[1])," *
   1487                           "$(sim.grains[i1[i]].lin_pos[2]) to " *
   1488                           "$(sim.grains[i2[i]].lin_pos[1])," *
   1489                           "$(sim.grains[i2[i]].lin_pos[2]) ")
   1490                     if tensile_stress[i] > 0
   1491                         write(f, "nohead ")
   1492                     else
   1493                         write(f, "heads ")
   1494                     end
   1495                     write(f, "lw $(abs(tensile_stress[i])/
   1496                                    max_tensile_stress*max_line_width) ")
   1497                     write(f, "lc palette cb $(tensile_stress[i])\n")
   1498                 end
   1499             end
   1500         end
   1501 
   1502         #write(f,"""plot "$(datafile)" with circles lt 1 lc rgb "black" t "Particle"
   1503         write(f,"""plot "$(datafile)" with circles fill solid fillcolor palette lw 0 title "Particle"
   1504               """)
   1505     end
   1506 
   1507     try
   1508         run(`gnuplot $gnuplotscript`)
   1509     catch return_signal
   1510         if isa(return_signal, Base.IOError)
   1511             error("Could not launch external gnuplot process")
   1512         end
   1513     end
   1514 
   1515     if verbose
   1516         @info filename
   1517     end
   1518 
   1519     if show_figure
   1520         if Sys.isapple()
   1521             run(`open $(filename)`)
   1522         elseif Sys.islinux()
   1523             run(`xdg-open $(filename)`)
   1524         end
   1525     end
   1526 end