tio.jl - Granular.jl - Julia package for granular dynamics simulation
HTML git clone git://src.adamsgaard.dk/Granular.jl
DIR Log
DIR Files
DIR Refs
DIR README
DIR 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 `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