diff --git a/src/Models/Material/Material_Factory.jl b/src/Models/Material/Material_Factory.jl
index 2e04b096..ddb19b8e 100644
--- a/src/Models/Material/Material_Factory.jl
+++ b/src/Models/Material/Material_Factory.jl
@@ -7,7 +7,7 @@ include("../../Core/Module_inclusion/set_Modules.jl")
 include("material_basis.jl")
 using LinearAlgebra
 using .Material_Basis:
-    get_all_elastic_moduli, distribute_forces, check_symmetry, get_all_elastic_moduli
+    get_all_elastic_moduli, distribute_forces!, check_symmetry, get_all_elastic_moduli
 using .Set_modules
 using TimerOutputs
 using StaticArrays
@@ -246,7 +246,8 @@ function distribute_force_densities(
     if !isnothing(nlist_filtered_ids)
         bond_norm = datamanager.get_field("Bond Norm")
         displacements = datamanager.get_field("Displacements", "NP1")
-        force_densities = distribute_forces(
+        force_densities = distribute_forces!(
+            force_densities,
             nodes,
             nlist,
             nlist_filtered_ids,
@@ -255,10 +256,9 @@ function distribute_force_densities(
             bond_damage,
             displacements,
             bond_norm,
-            force_densities,
         )
     else
-        distribute_forces(nodes, nlist, bond_force, volume, bond_damage, force_densities)
+        distribute_forces!(force_densities, nodes, nlist, bond_force, volume, bond_damage)
     end
 
 end
diff --git a/src/Models/Material/material_basis.jl b/src/Models/Material/material_basis.jl
index 3e1b3a88..53c3b871 100644
--- a/src/Models/Material/material_basis.jl
+++ b/src/Models/Material/material_basis.jl
@@ -11,7 +11,7 @@ using .Helpers: get_MMatrix, determinant, invert, smat
 export get_value
 export get_all_elastic_moduli
 export get_Hooke_matrix
-export distribute_forces
+export distribute_forces!
 export flaw_function
 export matrix_to_voigt
 export voigt_to_matrix
@@ -275,7 +275,7 @@ function get_Hooke_matrix(parameter::Dict, symmetry::String, dof::Int64, ID::Int
 end
 
 """
-    distribute_forces(nodes::Union{SubArray,Vector{Int64}}, nlist::Vector{Vector{Int64}}, nlist_filtered_ids::Vector{Vector{Int64}}, bond_force::Vector{Matrix{Float64}}, volume::Vector{Float64}, bond_damage::Vector{Vector{Float64}}, displacements::Matrix{Float64}, bond_norm::Vector{Matrix{Float64}}, force_densities::Matrix{Float64})
+    distribute_forces!(nodes::Union{SubArray,Vector{Int64}}, nlist::Vector{Vector{Int64}}, nlist_filtered_ids::Vector{Vector{Int64}}, bond_force::Vector{Matrix{Float64}}, volume::Vector{Float64}, bond_damage::Vector{Vector{Float64}}, displacements::Matrix{Float64}, bond_norm::Vector{Matrix{Float64}}, force_densities::Matrix{Float64})
 
 Distribute the forces on the nodes
 
@@ -292,7 +292,8 @@ Distribute the forces on the nodes
 # Returns
 - `force_densities::Matrix{Float64}`: The force densities.
 """
-function distribute_forces(
+function distribute_forces!(
+    force_densities::Matrix{Float64},
     nodes::Union{SubArray,Vector{Int64}},
     nlist::Vector{Vector{Int64}},
     nlist_filtered_ids::Vector{Vector{Int64}},
@@ -301,7 +302,6 @@ function distribute_forces(
     bond_damage::Vector{Vector{Float64}},
     displacements::Matrix{Float64},
     bond_norm::Vector{Vector{Vector{Float64}}},
-    force_densities::Matrix{Float64},
 )
 
     @inbounds @fastmath for iID in nodes
@@ -339,7 +339,7 @@ function distribute_forces(
 end
 
 """
-    distribute_forces(nodes::Union{SubArray,Vector{Int64}}, nlist::Vector{Vector{Int64}}, bond_force::Vector{Matrix{Float64}}, volume::Vector{Float64}, bond_damage::Vector{Vector{Float64}}, force_densities::Matrix{Float64})
+    distribute_forces!(nodes::Union{SubArray,Vector{Int64}}, nlist::Vector{Vector{Int64}}, bond_force::Vector{Matrix{Float64}}, volume::Vector{Float64}, bond_damage::Vector{Vector{Float64}}, force_densities::Matrix{Float64})
 
 Distribute the forces on the nodes
 
@@ -353,13 +353,13 @@ Distribute the forces on the nodes
 # Returns
 - `force_densities::Matrix{Float64}`: The force densities.
 """
-function distribute_forces(
+function distribute_forces!(
+    force_densities::Matrix{Float64},
     nodes::Union{SubArray,Vector{Int64}},
     nlist::Vector{Vector{Int64}},
     bond_force::Vector{Vector{Vector{Float64}}},
     volume::Vector{Float64},
     bond_damage::Vector{Vector{Float64}},
-    force_densities::Matrix{Float64},
 )
     @inbounds @fastmath for iID in nodes
         @views @inbounds @fastmath for jID ∈ axes(nlist[iID], 1)
@@ -375,7 +375,6 @@ function distribute_forces(
             end
         end
     end
-    return force_densities
 end
 
 
diff --git a/src/Support/helpers.jl b/src/Support/helpers.jl
index 54035b24..b73a66d5 100644
--- a/src/Support/helpers.jl
+++ b/src/Support/helpers.jl
@@ -360,14 +360,14 @@ Checks if the sum of the array is finite. If not, an error is raised.
 - `Bool`: `true` if the sum of the array is finite, `false` otherwise.
 """
 function check_inf_or_nan(array, msg)
-    if !isfinite(sum(array))
-        @error "Field ''$msg'' is infinite."
-        return true
-    end
     if isnan(sum(array))
         @error "Field ''$msg'' has NaN elements."
         return true
     end
+    if !isfinite(sum(array))
+        @error "Field ''$msg'' is infinite."
+        return true
+    end
     return false
 end
 
diff --git a/test/fullscale_tests/test_DCB/Reference/DCBmodel_PD_solid.e b/test/fullscale_tests/test_DCB/Reference/DCBmodel_PD_solid.e
index ad5fccf5..8f1e7088 100644
Binary files a/test/fullscale_tests/test_DCB/Reference/DCBmodel_PD_solid.e and b/test/fullscale_tests/test_DCB/Reference/DCBmodel_PD_solid.e differ
diff --git a/test/unit_tests/Models/Material/ut_material_basis.jl b/test/unit_tests/Models/Material/ut_material_basis.jl
index 96276705..e60aaa30 100644
--- a/test/unit_tests/Models/Material/ut_material_basis.jl
+++ b/test/unit_tests/Models/Material/ut_material_basis.jl
@@ -8,7 +8,7 @@ using .Material_Basis:
     get_value,
     get_all_elastic_moduli,
     get_Hooke_matrix,
-    distribute_forces,
+    distribute_forces!,
     flaw_function,
     matrix_to_voigt,
     voigt_to_matrix,
@@ -72,9 +72,16 @@ end
             bond_damage[iID] .* mapreduce(permutedims, vcat, bond_force[iID]) .* volume[iID]
     end
 
-    result =
-        distribute_forces(nodes, nlist, bond_force, volume, bond_damage, force_densities)
-    @test result ≈ expected_force_densities
+    distribute_forces!(
+        force_densities,
+        nodes,
+        nlist,
+        bond_force,
+        volume,
+        bond_damage,
+        force_densities,
+    )
+    @test force_densities ≈ expected_force_densities
 end
 @testset "ut_flaw_function" begin