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