diff --git a/examples/converge_diverge_flow.jl b/examples/converge_diverge_flow.jl
index 5366788..00ebab9 100644
--- a/examples/converge_diverge_flow.jl
+++ b/examples/converge_diverge_flow.jl
@@ -1,5 +1,4 @@
-using JLD2, Random, SplitApplyCombine, Statistics, StructArrays, Subzero
-import LibGEOS as LG
+using JLD2, Random, Statistics, Subzero
 
 # User Inputs
 const FT = Float64
@@ -72,6 +71,9 @@ simulation = Simulation(
 # Run simulation
 run!(simulation)
 
-Subzero.create_sim_gif("output/converge_diverge/floes.jld2", 
-                       "output/converge_diverge/initial_state.jld2",
-                       "output/converge_diverge/converge_diverge.gif")
\ No newline at end of file
+plot_sim(
+    "output/converge_diverge/floes.jld2",
+    "output/converge_diverge/initial_state.jld2",
+    20,
+    "output/converge_diverge/converge_diverge.mp4",
+)
\ No newline at end of file
diff --git a/examples/shear_flow.jl b/examples/shear_flow.jl
index 8730a6f..bf0a6b0 100644
--- a/examples/shear_flow.jl
+++ b/examples/shear_flow.jl
@@ -1,5 +1,4 @@
-using JLD2, Random, Statistics, StructArrays, Subzero
-import LibGEOS as LG
+using JLD2, Random, Statistics, Subzero
 
 # User Inputs
 const FT = Float64
@@ -79,8 +78,9 @@ simulation = Simulation(
 )
 run_time!(simulation)
 
-#ProfileView.@profview run!(simulation)
- 
-Subzero.create_sim_gif("output/shear_flow/floes.jld2", 
-                       "output/shear_flow/initial_state.jld2",
-                       "output/shear_flow/shear_flow.gif")
\ No newline at end of file
+plot_sim(
+    "output/shear_flow/floes.jld2",
+    "output/shear_flow/initial_state.jld2",
+    20,
+    "output/shear_flow/shear_flow.mp4",
+)
\ No newline at end of file
diff --git a/examples/simple_strait.jl b/examples/simple_strait.jl
index 29cfe2b..e2af5f9 100644
--- a/examples/simple_strait.jl
+++ b/examples/simple_strait.jl
@@ -1,13 +1,5 @@
-using JLD2, Random, SplitApplyCombine, Statistics, StructArrays, Subzero
-import LibGEOS as LG
+using JLD2, Random, Statistics, Subzero
 
-
-Subzero.plot_sim(
-    "output/simple_strait/floes.jld2",
-    "output/simple_strait/initial_state.jld2",
-    20,
-    "output/simple_strait/simple_strait.mp4",
-)
 # User Inputs
 const FT = Float64
 const Lx = 1e5
@@ -52,7 +44,7 @@ coupling_settings = CouplingSettings(
 # Floe creation
 floe_arr = initialize_floe_field(
     FT,
-    200,
+    75,
     [0.7],
     domain,
     hmean,
@@ -72,7 +64,7 @@ fracture_settings = FractureSettings(
         criteria = HiblerYieldCurve(floe_arr),
         Δt = 75,
         npieces = 3,
-        nhistory = 1000,
+        nhistory = 100,
         deform_on = false,
 )
 
@@ -88,7 +80,7 @@ simulation = Simulation(
     model = model,
     consts = consts,
     Δt = Δt,
-    nΔt = 20000,
+    nΔt = 5000,
     verbose = true,
     coupling_settings = coupling_settings,
     fracture_settings = fracture_settings,
@@ -97,9 +89,8 @@ simulation = Simulation(
     
 # Run simulation
 run_time!(simulation)
-#ProfileView.@profview run!(simulation)
 
-Subzero.plot_sim(
+plot_sim(
     "output/simple_strait/floes.jld2",
     "output/simple_strait/initial_state.jld2",
     Δt,
diff --git a/examples/test_run.jl b/examples/test_run.jl
index a40cca7..f97ff30 100644
--- a/examples/test_run.jl
+++ b/examples/test_run.jl
@@ -1,6 +1,69 @@
-using JLD2, Random, Statistics, Subzero, BenchmarkTools, StructArrays
+using JLD2, Random, Statistics, Subzero, BenchmarkTools, StructArrays, SplitApplyCombine
 import LibGEOS as LG
 
+
+Δt = 10
+Lx = 1e5
+Ly = 1e5
+collision_settings = CollisionSettings()
+grid = RegRectilinearGrid(
+    (-Lx, Lx),
+    (-Lx, Lx),
+    1e4,
+    1e4,
+)
+double_periodic_domain = Domain(
+    PeriodicBoundary(North, grid),
+    PeriodicBoundary(South, grid),
+    PeriodicBoundary(East, grid),
+    PeriodicBoundary(West, grid),
+)
+# Parent-parent collison (parents are touching)
+coords1 = splitdims(vcat([5*Lx/8 5*Lx/8 3*Lx/4 3*Lx/4].+1000, [3*Ly/4 5*Ly/4 5*Ly/4 3*Ly/4]))
+coords2 = splitdims(vcat(-[5*Lx/4 5*Lx/4 3*Lx/4-1000 3*Lx/4-1000], -[7*Lx/8 3*Lx/4-1000 3*Lx/4-1000 7*Lx/8]))
+
+floe_arr = StructArray(Floe([c], 0.5, 0.0) for c in [coords1, coords2])
+for i in eachindex(floe_arr)
+    floe_arr.id[i] = Float64(i)
+end
+trans_arr = StructArray([
+    Floe(
+        Subzero.translate([coords1],
+        0.0, -2Ly),
+        0.5,
+        0.0,
+    ),
+    Floe(
+        Subzero.translate([coords2], 2Lx, 0.0),
+        0.5,
+        0.0,
+    ),
+])
+for i in eachindex(trans_arr)
+    trans_arr.id[i] = i
+end
+spinlock = Threads.SpinLock()
+Subzero.timestep_collisions!(
+    trans_arr,
+    2,
+    double_periodic_domain,
+    Subzero.Constants(),
+    Δt,
+    collision_settings,
+    spinlock,
+)
+add_ghosts!(floe_arr, double_periodic_domain)
+Subzero.timestep_collisions!(
+    floe_arr,
+    2,
+    double_periodic_domain,
+    Subzero.Constants(),
+    Δt,
+    collision_settings,
+    spinlock,
+)
+
+
 # User Inputs
 const FT = Float64
 const Lx = 1e5
diff --git a/examples/uniform_flow.jl b/examples/uniform_flow.jl
index 85d19a4..10c74fa 100644
--- a/examples/uniform_flow.jl
+++ b/examples/uniform_flow.jl
@@ -64,7 +64,10 @@ simulation = Simulation(
     coupling_settings = CouplingSettings(two_way_coupling_on = true),
 )
 run_time!(simulation)
- 
-Subzero.create_sim_gif("output/uniform_flow/floes.jld2", 
-                       "output/uniform_flow/initial_state.jld2",
-                       "output/uniform_flow/uniform_flow.gif")
\ No newline at end of file
+
+Subzero.plot_sim(
+    "output/uniform_flow/floes.jld2",
+    "output/uniform_flow/initial_state.jld2",
+    20,
+    "output/uniform_flow/uniform_flow.mp4",
+)
\ No newline at end of file
diff --git a/src/Subzero.jl b/src/Subzero.jl
index 8186d09..1e36009 100644
--- a/src/Subzero.jl
+++ b/src/Subzero.jl
@@ -62,11 +62,12 @@ export
     IceStressCell,
     CellFloes,
     MonteCarloPointsGenerator,
-    SubGridPointsGenerator
+    SubGridPointsGenerator,
+    plot_sim
 
 import Base.@kwdef # this is being exported as of version 1.9
 import LibGEOS as LG
-using DataStructures, Dates, GeometryBasics, Interpolations, JLD2,
+using CairoMakie, DataStructures, Dates, GeometryBasics, Interpolations, JLD2,
     LinearAlgebra, Logging, Makie, Measures, NCDatasets, NetCDF,
     PolygonInbounds, Printf, Random, SplitApplyCombine, StaticArrays,
     Statistics, StructArrays, VoronoiCells
diff --git a/src/physical_processes/collisions.jl b/src/physical_processes/collisions.jl
index c541543..d9cb63e 100644
--- a/src/physical_processes/collisions.jl
+++ b/src/physical_processes/collisions.jl
@@ -296,6 +296,32 @@ function calc_friction_forces(
     return force
 end
 
+function add_interactions!(np, ifloe, idx::FT, forces, points, overs) where FT
+    inter_spots = size(ifloe.interactions, 1) - ifloe.num_inters
+    @views for i in 1:np
+        if forces[i, 1] != 0 || forces[i, 2] != 0
+            ifloe.num_inters += 1
+            if inter_spots < 1
+                ifloe.interactions = vcat(
+                    ifloe.interactions,
+                    zeros(FT, np - i + 1 , 7)
+                )
+                inter_spots += np - i
+            end
+            ifloe.interactions[ifloe.num_inters, floeidx] = idx
+            ifloe.interactions[ifloe.num_inters, xforce:yforce] .=
+                forces[i, :]
+            ifloe.interactions[ifloe.num_inters, xpoint:ypoint] .=
+                points[i, :]
+            ifloe.interactions[ifloe.num_inters, torque] = FT(0)
+            ifloe.interactions[ifloe.num_inters, overlap] = overs[i]
+            ifloe.overarea += overs[i]
+            inter_spots -= 1
+        end
+    end
+    return
+end
+
 """
     floe_floe_interaction!(
         ifloe,
@@ -337,7 +363,6 @@ function floe_floe_interaction!(
     i,
     jfloe,
     j,
-    nfloes,
     consts,
     Δt,
     max_overlap::FT,
@@ -351,7 +376,6 @@ function floe_floe_interaction!(
         total_area += a
     end
     if total_area > 0
-        # Floes overlap too much - remove floe or transfer floe mass
         # Floes overlap too much - mark floes to be fused
         if max(total_area/ifloe.area, total_area/jfloe.area) > max_overlap
             ifloe.status.tag = fuse
@@ -390,16 +414,7 @@ function floe_floe_interaction!(
                 )
                 # Calculate total forces and update ifloe's interactions
                 forces = normal_forces .+ friction_forces
-                if sum(abs.(forces)) != 0
-                    new_interactions = [
-                        fill(j, np) forces fpoints zeros(FT, np) overlaps
-                    ]
-                    ifloe.interactions = vcat(
-                        ifloe.interactions,
-                        new_interactions
-                    )
-                    ifloe.overarea += sum(overlaps)
-                end
+                add_interactions!(np, ifloe, j, forces, fpoints, overlaps)
             end
         end
     end
@@ -650,11 +665,7 @@ function floe_domain_element_interaction!(
                 )
                 # Calculate total forces and update ifloe's interactions
                 forces = normal_forces .+ friction_forces
-                if sum(abs.(forces)) != 0
-                    floe.interactions = [floe.interactions;
-                        fill(Inf, np) forces fpoints zeros(np) overlaps]
-                    floe.overarea += sum(overlaps)
-                end
+                add_interactions!(np, floe, FT(Inf), forces, fpoints, overlaps)
             end
         end
     end
@@ -809,7 +820,6 @@ function floe_domain_interaction!(
             )
         end
     end
-
     return
 end
 
@@ -826,15 +836,13 @@ function calc_torque!(
     floe::Union{LazyRow{<:Floe{FT}}, Floe{FT}},
 ) where {FT<:AbstractFloat}
     inters = floe.interactions
-    if !isempty(inters)
-        dir = [inters[:, xpoint] .- floe.centroid[1] inters[:, ypoint] .- floe.centroid[2] zeros(FT, size(inters, 1))]
-        frc = [inters[:, xforce] inters[:, yforce] zeros(FT, size(inters, 1))]
-        for i in axes(dir, 1)
-            idir = vec(dir[i, :])
-            ifrc = vec(frc[i, :])
-            itorque = cross(idir, ifrc)
-            floe.interactions[i, torque] = itorque[3]
-        end
+    for i in 1:floe.num_inters
+        xp = inters[i, xpoint] .- floe.centroid[1]
+        yp = inters[i, ypoint] .- floe.centroid[2]
+        xf = inters[i, xforce]
+        yf = inters[i, yforce]
+        # cross product of direction and force where third element of both is 0
+        floe.interactions[i, torque] = xp * yf - yp * xf
     end
     return
 end
@@ -901,7 +909,7 @@ function timestep_collisions!(
         # reset collision values
         fill!(floes.collision_force[i], FT(0))
         floes.collision_trq[i] = FT(0.0)
-        floes.interactions[i] = zeros(FT, 0, 7)
+        floes.num_inters[i] = 0
         for j in i+1:length(floes)
             id_pair, ghost_id_pair =
                 if floes.id[i] > floes.id[j]
@@ -933,7 +941,6 @@ function timestep_collisions!(
                         i,
                         LazyRow(floes, j),
                         j,
-                        n_init_floes,
                         consts,
                         Δt,
                         collision_settings.floe_floe_max_overlap,
@@ -951,7 +958,7 @@ function timestep_collisions!(
     end
     # Move compression boundaries if they exist
     update_boundaries!(domain, Δt)
-    # Update floes not directly calculated above where i>j - can't be parallelized
+    # Update floes not directly calculated above where i>j - can't be parallel
     for i in eachindex(floes)
         # Update fuse information
         if floes.status[i].tag == fuse
@@ -961,39 +968,60 @@ function timestep_collisions!(
             end
         end
         # Update interaction information
-        ij_inters = floes.interactions[i]
-        if !isempty(ij_inters)
+        i_inters = floes.interactions[i]
+        if floes.num_inters[i] > 0
             # Loop over each interaction with floe i
-            for inter_idx in axes(ij_inters, 1)
-                j = ij_inters[inter_idx, floeidx]  # Index of floe to update
+            for inter_idx in 1:floes.num_inters[i]
+                j = i_inters[inter_idx, floeidx]  # Index of floe to update
+                # not boundaries and hasn't been updated yet
                 if j <= length(floes) && j > i
                     jidx = Int(j)
-                    floes.interactions[jidx] = vcat(
-                        floes.interactions[jidx],
-                        ij_inters[inter_idx:inter_idx, :]
+                    # add matching interaction (j with i)
+                    add_interactions!(1, LazyRow(floes, jidx), FT(i),
+                        i_inters[inter_idx:inter_idx, xforce:yforce],
+                        i_inters[inter_idx:inter_idx, xpoint:ypoint],
+                        i_inters[inter_idx:inter_idx, overlap:overlap]
                     )
-                    floes.interactions[jidx][end, floeidx] = i
-                    floes.interactions[jidx][end, xforce] *= -1
-                    floes.interactions[jidx][end, yforce] *= -1
-                    # floes.interactions[jidx] = [floes.interactions[jidx]; i -ij_inters[inter_idx, xforce] -ij_inters[inter_idx, yforce] #=
-                    #                          =# ij_inters[inter_idx, xpoint] ij_inters[inter_idx, ypoint] FT(0.0) ij_inters[inter_idx, overlap]]
-                    floes.overarea[jidx] += ij_inters[inter_idx, overlap]
+                    # equal and opposite forces
+                    j_np = floes.num_inters[jidx]
+                    floes.interactions[jidx][j_np, xforce:yforce] *= -1
                 end
             end
         end
     end
     # Calculate total on parents by summing interactions on parent and children
-    for i in range(1, n_init_floes)
-        for g in floes.ghosts[i]
-            g_inters = floes.interactions[g]
-            g_inters[:, xpoint] .-= (floes.centroid[g][1] - floes.centroid[i][1])
-            g_inters[:, ypoint] .-= (floes.centroid[g][2] - floes.centroid[i][2])
-            floes.interactions[i] = [floes.interactions[i]; g_inters]
+    for i in 1:n_init_floes
+        if !isempty(floes.ghosts[i])
+            for g in floes.ghosts[i]
+                gnp = floes.num_inters[g]
+                # shift interaction points to parent locations
+                floes.interactions[g][1:gnp, xpoint] .-=
+                    floes.centroid[g][1] - floes.centroid[i][1]
+                floes.interactions[g][1:gnp, ypoint] .-=
+                    floes.centroid[g][2] - floes.centroid[i][2]
+                # add interactions
+                add_interactions!(gnp, LazyRow(floes, i), i,
+                    floes.interactions[g][1:gnp, xforce:yforce],
+                    floes.interactions[g][1:gnp, xpoint:ypoint],
+                    floes.interactions[g][1:gnp, overlap]
+                )
+                inp = floes.num_inters[i]
+                # set interaction indices to correct values
+                floes.interactions[i][(inp - gnp + 1):inp, floeidx] .= 
+                    floes.interactions[g][1:gnp, floeidx]
+            end
         end
+        # calculate interactions torque and total forces / torque
         calc_torque!(LazyRow(floes, i))
-        floes.collision_force[i][1] += sum(@view floes.interactions[i][:, xforce])
-        floes.collision_force[i][2] += sum(@view floes.interactions[i][:, yforce])
-        floes.collision_trq[i] += sum(@view floes.interactions[i][:, torque])
+        floes.collision_force[i][1] += sum(
+            @view floes.interactions[i][1:floes.num_inters[i], xforce]
+        )
+        floes.collision_force[i][2] += sum(
+            @view floes.interactions[i][1:floes.num_inters[i], yforce]
+        )
+        floes.collision_trq[i] += sum(
+            @view floes.interactions[i][1:floes.num_inters[i], torque]
+        )
     end
     return
 end
diff --git a/src/physical_processes/fractures.jl b/src/physical_processes/fractures.jl
index 222a5ec..a29867c 100644
--- a/src/physical_processes/fractures.jl
+++ b/src/physical_processes/fractures.jl
@@ -524,7 +524,8 @@ function fracture_floes!(
     Threads.@threads for i in 1:nfloes2frac
         # Deform floe around largest impact site
         if fracture_settings.deform_on
-            inters = floes.interactions[frac_idx[i]]
+            fidx = frac_idx[i]
+            inters = floes.interactions[fidx][1:floes.num_inters[fidx], :]
             inters = inters[.!(isinf.(inters[:, floeidx])), :]
             if !isempty(inters)
                 _, max_inters_idx = findmax(inters[:, overlap])
diff --git a/src/physical_processes/update_floe.jl b/src/physical_processes/update_floe.jl
index 160b16d..4d51450 100644
--- a/src/physical_processes/update_floe.jl
+++ b/src/physical_processes/update_floe.jl
@@ -239,17 +239,21 @@ Inputs:
 Outputs:
     Caculates stress on floe at current timestep from interactions
 """
-function calc_stress!(floe)
+function calc_stress!(floe::Union{LazyRow{Floe{FT}}, Floe{FT}}) where {FT}
     # Stress calcultions
     xi, yi = floe.centroid
     inters = floe.interactions
     # Calculates timestep stress
-    stress = fill(1/(2* floe.area * floe.height), 2, 2)
-    stress[1, 1] *= 2sum((inters[:, xpoint] .- xi) .* inters[:, xforce])
-    stress[1, 2] *= sum((inters[:, ypoint] .- yi) .* inters[:, xforce]) + 
-        sum(inters[:, yforce] .* (inters[:, xpoint] .- xi))
+    stress = fill(FT(0), 2, 2)
+    for i in 1:floe.num_inters
+        stress[1, 1] += (inters[i, xpoint] - xi) * inters[i, xforce]
+        stress[1, 2] += (inters[i, ypoint] - yi) * inters[i, xforce] +
+            (inters[i, xpoint] - xi) * inters[i, yforce]
+        stress[2, 2] += (inters[i, ypoint] - yi) * inters[i, yforce]
+    end
+    stress[1, 2] *= FT(0.5)
     stress[2, 1] = stress[1, 2]
-    stress[2, 2] *= 2sum((inters[:, ypoint] .- yi) .* inters[:, yforce])
+    stress .*= 1/(floe.area * floe.height)
     # Add timestep stress to stress history
     push!(floe.stress_history, stress)
     # Average stress history to find floe's average stress
@@ -266,35 +270,31 @@ Inputs:
 Outputs:
     strain      <Matrix{AbstractFloat}> 2x2 matrix for floe strain 
 """
-function calc_strain!(floe)
+function calc_strain!(floe::Union{LazyRow{Floe{FT}}, Floe{FT}}) where {FT}
     # coordinates of floe centered at centroid
     translate!(floe.coords, -floe.centroid[1], -floe.centroid[2])
-    xcoords, ycoords = separate_xy(floe.coords)
-    translate!(floe.coords, floe.centroid[1], floe.centroid[2])
-    # Find distance between each vertex
-    if xcoords[1] != xcoords[end] && ycoords[1] != ycoords[end]
-        push!(xcoords, xcoords[1])
-        push!(ycoords, ycoords[1])
+    fill!(floe.strain, FT(0))
+    for i in 1:(length(floe.coords[1]) - 1)
+        xdiff = floe.coords[1][i + 1][1] - floe.coords[1][i][1]
+        ydiff = floe.coords[1][i + 1][2] - floe.coords[1][i][2]
+        rad1 = sqrt(floe.coords[1][i][1]^2 + floe.coords[1][i][2]^2)
+        θ1 = atan(floe.coords[1][i][2], floe.coords[1][i][1])
+        rad2 = sqrt(floe.coords[1][i + 1][1]^2 + floe.coords[1][i + 1][2]^2)
+        θ2 = atan(floe.coords[1][i + 1][2], floe.coords[1][i + 1][1])
+        u1 = floe.u - floe.ξ * rad1 * sin(θ1)
+        u2 = floe.u - floe.ξ * rad2 * sin(θ2)
+        v1 = floe.u + floe.ξ * rad1 * cos(θ1)
+        v2 = floe.u + floe.ξ * rad2 * cos(θ2)
+        udiff = u2 - u1
+        vdiff = v2 - v1
+        floe.strain[1, 1] += udiff * ydiff
+        floe.strain[1, 2] += udiff * xdiff + vdiff * ydiff
+        floe.strain[2, 2] += vdiff * xdiff
     end
-    xcoords_diff = diff(xcoords)
-    ycoords_diff = diff(ycoords)
-    # u and v velocities of floes at each vertex
-    ucoords = fill(floe.u, size(xcoords))
-    vcoords = fill(floe.v, size(xcoords))
-    for i in eachindex(ucoords)
-        rad = sqrt(xcoords[i]^2 + ycoords[i]^2)
-        θ = atan(ycoords[i], xcoords[i])
-        ucoords[i] -= floe.ξ * rad * sin(θ)
-        vcoords[i] += floe.ξ * rad * cos(θ)
-    end
-    ucoords_diff = diff(ucoords)
-    vcoords_diff = diff(vcoords)
-    fill!(floe.strain, 1/(2 * floe.area))
-    floe.strain[1, 1] *= sum(ucoords_diff .* ycoords_diff) # dudx
-    floe.strain[1, 2] *= 0.5(sum(ucoords_diff .* xcoords_diff) +
-        sum(vcoords_diff .* ycoords_diff)) # dudy + dvdx
+    floe.strain[1, 2] *= FT(0.5)
     floe.strain[2, 1] = floe.strain[1, 2]
-    floe.strain[2, 2] *= sum(vcoords_diff .* xcoords_diff) # dvdy
+    floe.strain ./= 2floe.area
+    translate!(floe.coords, floe.centroid[1], floe.centroid[2])
     return
 end
 
@@ -322,7 +322,7 @@ function timestep_floe_properties!(
         cforce = floes.collision_force[i]
         ctrq = floes.collision_trq[i]
         # Update stress
-        if !isempty(floes.interactions[i])
+        if floes.num_inters[i] > 0
             calc_stress!(LazyRow(floes, i))
         end
         # Ensure no extreem values due to model instability
diff --git a/src/simulation_components/floe.jl b/src/simulation_components/floe.jl
index b0f0cd6..84e853e 100644
--- a/src/simulation_components/floe.jl
+++ b/src/simulation_components/floe.jl
@@ -3,7 +3,7 @@ Structs and functions used to define floes and floe fields within Subzero
 """
 
 """
-Enum to index into floe interactions field with more intuituve names
+Enum for differnt floe status
 """
 @enum StatusTag begin
     active = 1
@@ -136,6 +136,7 @@ Singular sea ice floe with fields describing current state.
     collision_force::Matrix{FT} = zeros(1, 2)
     collision_trq::FT = 0.0
     interactions::Matrix{FT} = zeros(0, 7)
+    num_inters::Int = 0
     stress::Matrix{FT} = zeros(2, 2)
     stress_history::StressCircularBuffer{FT} = StressCircularBuffer(1000)
     strain::Matrix{FT} = zeros(2, 2)
diff --git a/src/simulation_components/simulation.jl b/src/simulation_components/simulation.jl
index 325c82d..adce441 100644
--- a/src/simulation_components/simulation.jl
+++ b/src/simulation_components/simulation.jl
@@ -141,8 +141,6 @@ function timestep_sim!(sim, tstep)
             sim.Δt,
             sim.simp_settings.max_floe_height,
         )
-
-        # TODO: Remove parent ids ?
         # Fracture floes
         if sim.fracture_settings.fractures_on && mod(tstep, sim.fracture_settings.Δt) == 0
             max_floe_id =
diff --git a/test/inputs/stress_strain.jld2 b/test/inputs/stress_strain.jld2
new file mode 100644
index 0000000..72d2eb7
Binary files /dev/null and b/test/inputs/stress_strain.jld2 differ
diff --git a/test/inputs/test_values.jld2 b/test/inputs/test_values.jld2
deleted file mode 100644
index 8017e9e..0000000
Binary files a/test/inputs/test_values.jld2 and /dev/null differ
diff --git a/test/qualitative_behavior.jl b/test/qualitative_behavior.jl
index ee09139..36845be 100644
--- a/test/qualitative_behavior.jl
+++ b/test/qualitative_behavior.jl
@@ -12,6 +12,7 @@ const Δgrid = 10000
 const hmean = 0.25
 const Δh = 0.0
 const nΔt = 4000
+const Δt = 10
 const newfloe_Δt = 500
 const coarse_nx = 10
 const coarse_ny = 10
@@ -126,7 +127,7 @@ writers1 = OutputWriters(
 simulation1 = Simulation(
     name = "sim1",
     model = model1,
-    Δt = 10,
+    Δt = Δt,
     nΔt = nΔt,
     collision_settings = collisions_off_settings,
     writers = writers1,
@@ -161,7 +162,7 @@ writers2 = OutputWriters(
 simulation2 = Simulation(
     name = "sim2",
     model = model2,
-    Δt = 10,
+    Δt = Δt,
     nΔt = nΔt,
     collision_settings = collisions_off_settings,
     writers = writers2,
@@ -198,7 +199,7 @@ writers3 = OutputWriters(
 simulation3 = Simulation(
     name = "sim3",
     model = model3,
-    Δt = 10,
+    Δt = Δt,
     nΔt = nΔt,
     writers = writers3,
     coupling_settings = CouplingSettings(coupling_on = false),
@@ -274,7 +275,7 @@ writers4 = OutputWriters(
 simulation4 = Simulation(
     name = "sim4",
     model = model4,
-    Δt = 10,
+    Δt = Δt,
     nΔt = nΔt,
     writers = writers4,
     coupling_settings = CouplingSettings(coupling_on = false),
@@ -323,7 +324,7 @@ writers5 = OutputWriters(
 simulation5 = Simulation(
     name = "sim5",
     model = model5,
-    Δt = 10,
+    Δt = Δt,
     nΔt = nΔt,
     writers = writers5,
     coupling_settings = CouplingSettings(coupling_on = false),
@@ -335,9 +336,10 @@ push!(sim_arr, simulation5)
 # Run the simulations
 for sim in sim_arr
     run!(sim)
-    Subzero.create_sim_gif(
+    plot_sim(
         joinpath("test/output", sim.name, "floes.jld2"),
         joinpath("test/output", sim.name, "initial_state.jld2"),
-        joinpath("test/output", sim.name, string(sim.name, ".gif")),
+        Δt,
+        joinpath("test/output", sim.name, string(sim.name, ".mp4")),
     )
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 45c3125..ba5aa0d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,16 +5,16 @@ import LibGEOS as LG
 using Test
 
 @testset "Subzero.jl" begin
-    # include("test_physical_processes/test_update_floe.jl")
+    include("test_physical_processes/test_update_floe.jl")
     include("test_physical_processes/test_collisions.jl")
-    # include("test_physical_processes/test_coupling.jl")
-    # include("test_physical_processes/test_process_settings.jl")
-    # include("test_physical_processes/test_fractures.jl")
-    # include("test_physical_processes/test_simplification.jl")
-    # include("test_floe.jl")
-    # include("test_floe_utils.jl")
-    # include("test_model.jl")
-    # include("test_output.jl")
-    # include("test_simulation.jl")
-    # include("test_conservation.jl")
+    include("test_physical_processes/test_coupling.jl")
+    include("test_physical_processes/test_process_settings.jl")
+    include("test_physical_processes/test_fractures.jl")
+    include("test_physical_processes/test_simplification.jl")
+    include("test_floe.jl")
+    include("test_floe_utils.jl")
+    include("test_model.jl")
+    include("test_output.jl")
+    include("test_simulation.jl")
+    include("test_conservation.jl")
 end
diff --git a/test/test_physical_processes/test_collisions.jl b/test/test_physical_processes/test_collisions.jl
index 6f2d623..bbb56fc 100644
--- a/test/test_physical_processes/test_collisions.jl
+++ b/test/test_physical_processes/test_collisions.jl
@@ -41,7 +41,7 @@
         cfloe.u = 0.3
         consts = Constants()
         # Triange tip intersected with a rectangle
-        Subzero.floe_floe_interaction!(tri, 1, rect, 2, 2, consts, Δt, max_overlap)
+        Subzero.floe_floe_interaction!(tri, 1, rect, 2, consts, Δt, max_overlap)
         @test isapprox(tri.interactions[1, xforce], -64613382.47, atol = 1e-2)
         @test isapprox(tri.interactions[1, yforce], -521498991.51, atol = 1e-2)
         @test isapprox(tri.interactions[1, xpoint], 10000.00, atol = 1e-2)
@@ -53,7 +53,7 @@
         @test isapprox(tri.interactions[1, torque], 1069710443203.99, atol = 1e-2)
         tri.interactions = zeros(0, 7)
         # Sideways C intersected with rectangle so there are two areas of overlap
-        Subzero.floe_floe_interaction!(cfloe, 1, rect, 2, 2, consts, Δt, max_overlap)
+        Subzero.floe_floe_interaction!(cfloe, 1, rect, 2, consts, Δt, max_overlap)
         @test isapprox(cfloe.interactions[1, xforce], -163013665.41, atol = 1e-2)
         @test isapprox(cfloe.interactions[2, xforce], -81506832.70, atol = 1e-2)
         @test isapprox(cfloe.interactions[1, yforce], 804819565.60, atol = 1e-2)
@@ -72,7 +72,7 @@
         # Floes overlapping more than 55%  - rectangle and shifted rectangle
         shift_rect = deepcopy(rect)
         shift_rect.coords = Subzero.translate(shift_rect.coords, 0.5e4, 0.0)
-        Subzero.floe_floe_interaction!(rect, 1, shift_rect, 2, 2, consts, Δt, max_overlap)
+        Subzero.floe_floe_interaction!(rect, 1, shift_rect, 2, consts, Δt, max_overlap)
         @test rect.status.tag == Subzero.fuse 
         @test rect.status.fuse_idx == [2]
         @test isempty(rect.interactions)
@@ -95,7 +95,6 @@
             1,
             small_rect,
             2,
-            2,
             consts,
             Δt,
             max_overlap,
@@ -106,7 +105,7 @@
         # Overlapping barely floes - such a small overlap that forces are not calculated
         shift_rect = deepcopy(rect)
         shift_rect.coords = Subzero.translate(shift_rect.coords, 1.9999999e4, 0.0)
-        Subzero.floe_floe_interaction!(shift_rect, 1, rect, 2, 2, consts, Δt, max_overlap)
+        Subzero.floe_floe_interaction!(shift_rect, 1, rect, 2, consts, Δt, max_overlap)
         @test isempty(shift_rect.interactions)
     end
 
@@ -252,6 +251,7 @@
         # Test floe overlapping >75% with collision boundary -> but max overlap is now 1
         Subzero.floe_domain_interaction!(efloe_large, domain, consts, Δt, 1.0)
         @test !isempty(efloe_large.interactions)
+        @test efloe_large.num_inters > 0
         # Test floe passing through open boundary is killed
         Subzero.floe_domain_interaction!(wfloe, domain, consts, Δt, max_overlap)
         @test wfloe.status.tag == Subzero.remove
diff --git a/test/test_physical_processes/test_fractures.jl b/test/test_physical_processes/test_fractures.jl
index 9ce6836..84aac2b 100644
--- a/test/test_physical_processes/test_fractures.jl
+++ b/test/test_physical_processes/test_fractures.jl
@@ -156,6 +156,7 @@
             -148920620521.112,
             6795329.38154967,
         ]')
+        frac_deform_floe.num_inters = 1
         frac_deform_floe.p_dudt = 0.11
         frac_floe.stress = frac_stress
         no_frac_small.stress = frac_stress
diff --git a/test/test_physical_processes/test_update_floe.jl b/test/test_physical_processes/test_update_floe.jl
index 2ea370e..1acd43b 100644
--- a/test/test_physical_processes/test_update_floe.jl
+++ b/test/test_physical_processes/test_update_floe.jl
@@ -21,15 +21,14 @@
 
         # Test Stress and Strain Calculations
         floe_dict = load(
-            "inputs/test_values.jld2"  # uses the first 2 element
+            "inputs/stress_strain.jld2"  # uses the first 2 element
         )
         stresses = [[-10.065, 36.171, 36.171, -117.458],
             [7.905, 21.913, 21.913, -422.242]]
         stress_histories = [[-4971.252, 17483.052, 17483.052, -57097.458],
             [4028.520, 9502.886, 9502.886, -205199.791]]
-        strains = [[-3.724, 0, 0, 0], [7.419, 0, 0,	-6.987]]
-        strain_multiplier = [1e28, 1e6]
-
+        strains = [[-0.0372, 0, 0, .9310], [7.419, 0, 0,	-6.987]]
+        strain_multiplier = [1e6, 1e6]
         for i in 1:2
             f = Floe(
                 floe_dict["coords"][i],
@@ -40,7 +39,8 @@
                 ξ = floe_dict["ξ"][i],
             )
             f.interactions = floe_dict["interactions"][i]
-            f.stress_history = floe_dict["stress_history"][i]
+            f.num_inters = size(f.interactions, 1)
+            push!(f.stress_history, floe_dict["last_stress"][i])
             stress = Subzero.calc_stress!(f)
             @test all(isapprox.(vec(f.stress), stresses[i], atol = 1e-3))
             @test all(isapprox.(
@@ -54,6 +54,7 @@
                 strains[i],
                 atol = 1e-3
             ))
+            @test f.coords == floe_dict["coords"][i]
         end
     end
     @testset "Replace floe" begin