GopherProxy

	tio.jl - Granular.jl - Julia package for granular dynamics simulation
	git clone git://src.adamsgaard.dk/Granular.jl
	Log
	Files
	Refs
	README
	LICENSE
	---
	tio.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 `fil…
	17 parameter is not specified, it will be saved to a subdirectory under the…
	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 fil…
	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 defa…
	68 value (`-1`), the function will try to read the latest file, determi…
	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.…
	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 fil…
	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.…
	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 fil…
	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…
	171 `loop=true`, as the script regenerates (and overwrites) all output …
	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.RFC1123F…
	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.060.024.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 ParaVi…
	293 can be visualized by applying a Glyph filter.
	294
	295 If the simulation contains an `Ocean` data structure, it's contents will…
	296 written to separate `.vtu` files. This can be disabled by setting the a…
	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 ParaVi…
	337 can be visualized by applying a Glyph filter. This function is called…
	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 …
	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 displacemen…
	373
	374 WriteVTK.vtk_point_data(vtkfile, ifarr.ang_pos, "Angular position [r…
	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 …
	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_displacemen…
	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 index…
	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=\"Float…
	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_…
	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 …
	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` i…
	743 structured grid (file type `.vts`). These files can be read by ParaView…
	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, verbos…
	822
	823 Create a `".py"` script for visualizing the simulation VTK files in Para…
	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 def…
	830 (blank) value, the script is named after the simulation id (`simulat…
	831 * `folder::String`: output directory, current directory the default.
	832 * `vtk_folder::String`: directory containing the VTK output files, by de…
	833 points to the full system path equivalent to `"./<simulation.id>/"`.
	834 * `save_animation::Bool`: make the generated script immediately save a r…
	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…
	942 imagegrainsDisplay.ScaleFactor = 6.050000000000001
	943 imagegrainsDisplay.SelectScaleArray = 'Angular acceleration [rad s^-2]'
	944 imagegrainsDisplay.GlyphType = 'Arrow'
	945 imagegrainsDisplay.GlyphTableIndexArray = 'Angular acceleration [rad s^-…
	946 imagegrainsDisplay.DataAxesGrid = 'GridAxesRepresentation'
	947 imagegrainsDisplay.PolarAxes = 'PolarAxesRepresentation'
	948 imagegrainsDisplay.ScalarOpacityUnitDistance = 64.20669746996803
	949 imagegrainsDisplay.GaussianRadius = 3.0250000000000004
	950 imagegrainsDisplay.SetScaleArray = ['POINTS', 'Atmosphere drag coefficie…
	951 imagegrainsDisplay.ScaleTransferFunction = 'PiecewiseFunction'
	952 imagegrainsDisplay.OpacityArray = ['POINTS', 'Atmosphere drag coefficien…
	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…
	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 (horizont…
	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…
	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…
	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 execut…
	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 defa…
	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, req…
	1103 `images=true`)
	1104 * `animation::Bool`: render animation as movie to disk (default: false).…
	1105 ffmpeg is available on the system, the `.avi` file is converted to `…
	1106 * `trim::Bool`: trim images in animated sequence (default: true)
	1107 * `reverse::Bool`: if `images=true` additionally render reverse-animated…
	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 availabl…
	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.…
	1205 plot is saved accoring to the simulation id, the optional `filename_post…
	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` r…
	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 (defaul…
	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…
	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)+b…
	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 (defaul…
	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 '#99…
	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…
	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" …
	1503 write(f,"""plot "$(datafile)" with circles fill solid fillcolor …
	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