Merge branch 'main' of https://github.com/PeriHub/PeriLab.jl

PeriHub · Nov 5, 2024 · 29aac55 · 29aac55
2 parents 405b9c4 + 02e49a9
commit 29aac55
Show file tree

Hide file tree

Showing 17 changed files with 219 additions and 137 deletions.
diff --git a/src/IO/mesh_data.jl b/src/IO/mesh_data.jl
@@ -365,8 +365,8 @@ function get_bond_geometry(datamanager::Module)
     undeformed_bond = datamanager.create_constant_bond_field("Bond Geometry", Float64, dof)
     undeformed_bond_length =
         datamanager.create_constant_bond_field("Bond Length", Float64, 1)
-    undeformed_bond, undeformed_bond_length = Geometry.bond_geometry(
-        Vector(1:nnodes),
+    Geometry.bond_geometry(
+        Vector{Int64}(1:nnodes),
         nlist,
         coor,
         undeformed_bond,

diff --git a/src/Models/Material/Material_Models/Correspondence/Bond_Associated_Correspondence.jl b/src/Models/Material/Material_Models/Correspondence/Bond_Associated_Correspondence.jl
@@ -16,7 +16,7 @@ using .Geometry:
     compute_bond_level_rotation_tensor,
     compute_bond_level_deformation_gradient
 include("../../../Pre_calculation/pre_bond_associated_correspondence.jl")
-using .Pre_Bond_Associated_Correspondence: compute_weighted_volume
+using .Pre_Bond_Associated_Correspondence: compute_weighted_volume!
 using TimerOutputs
 
 export init_model
@@ -69,8 +69,7 @@ function init_model(
     bond_damage = datamanager.get_bond_damage("NP1")
     weighted_volume = datamanager.create_constant_node_field("Weighted Volume", Float64, 1)
 
-    weighted_volume =
-        compute_weighted_volume(nodes, nlist, volume, bond_damage, omega, weighted_volume)
+    compute_weighted_volume!(weighted_volume, nodes, nlist, volume, bond_damage, omega)
 
     return datamanager
 end
@@ -175,7 +174,7 @@ function compute_model(
 
     ba_rotation_tensor = datamanager.get_field("Bond Rotation Tensor", "NP1")
 
-    strain_NP1, strain_increment = compute_bond_strain(
+    compute_bond_strain(
         nodes,
         nlist,
         ba_deformation_gradient,
@@ -336,20 +335,15 @@ function compute_bond_strain(
 )
 
     for iID in nodes
-        @views strain_NP1[iID][:, :, :] = compute_strain(
-            eachindex(nlist[iID]),
-            (deformation_gradient[iID][:, :, :]),
-            (strain_NP1[iID][:, :, :]),
-        )
-        @views matrix_diff!(
+        compute_strain(eachindex(nlist[iID]), deformation_gradient[iID], strain_NP1[iID])
+        matrix_diff!(
             strain_increment[iID],
             eachindex(nlist[iID]),
             strain_NP1[iID],
             strain_N[iID],
         )
 
     end
-    return strain_NP1, strain_increment
 end
 
 
@@ -437,10 +431,12 @@ function compute_bond_forces(
             end
 
 
-            bond_forces[iID][jID] .=
-                integral_nodal_stress[iID, :, :] * gradient_weights[iID][jID]
-            # mul!(bond_forces[iID][jID, :], integral_nodal_stress[iID, :, :], gradient_weights[iID][jID, :])
-            bond_forces[iID][jID] +=
+            # bond_forces[iID][jID, :] =
+            #     integral_nodal_stress[iID, :, :] * gradient_weights[iID][jID, :]
+            bf = @view bond_forces[iID][jID]
+            gw = @view gradient_weights[iID][jID]
+            mul!(bf, integral_nodal_stress[iID, :, :], gw)
+            @views bond_forces[iID][jID] +=
                 bond_damage[iID][jID] * omega[iID][jID] /
                 (weighted_volume[iID] * bond_length[iID][jID] * bond_length[iID][jID]) .*
                 bond_stress[iID][jID, :, :] * bond_geometry[iID][jID]

diff --git a/src/Models/Material/Material_Models/Correspondence/Correspondence.jl b/src/Models/Material/Material_Models/Correspondence/Correspondence.jl
@@ -190,8 +190,7 @@ function compute_model(
     stress_N = datamanager.get_field("Cauchy Stress", "N")
     stress_NP1 = datamanager.get_field("Cauchy Stress", "NP1")
     strain_increment = datamanager.get_field("Strain Increment")
-    strain_NP1 = compute_strain(nodes, deformation_gradient, strain_NP1)
-
+    compute_strain(nodes, deformation_gradient, strain_NP1)
     matrix_diff!(strain_increment, nodes, strain_NP1, strain_N)
 
     if rotation

diff --git a/src/Models/Material/Material_Models/Correspondence/Correspondence_Elastic.jl b/src/Models/Material/Material_Models/Correspondence/Correspondence_Elastic.jl
@@ -169,13 +169,10 @@ function compute_stresses_ba(
         @views hookeMatrix =
             get_Hooke_matrix(material_parameter, material_parameter["Symmetry"], dof, iID)
         @views for jID in eachindex(nlist[iID])
-            @views fast_mul!(
-                stress_NP1[iID][jID, :, :],
-                hookeMatrix,
-                strain_increment[iID][jID, :, :],
-                stress_N[iID][jID, :, :],
-                mapping,
-            )
+            @views sNP1 = stress_NP1[iID][jID, :, :]
+            @views sInc = strain_increment[iID][jID, :, :]
+            @views sN = stress_N[iID][jID, :, :]
+            fast_mul!(sNP1, hookeMatrix, sInc, sN, mapping)
         end
     end
     return stress_NP1, datamanager

diff --git a/src/Models/Pre_calculation/bond_deformation.jl b/src/Models/Pre_calculation/bond_deformation.jl
@@ -80,13 +80,7 @@ function compute(
     deformed_coor = datamanager.get_field("Deformed Coordinates", "NP1")
     deformed_bond = datamanager.get_field("Deformed Bond Geometry", "NP1")
     deformed_bond_length = datamanager.get_field("Deformed Bond Length", "NP1")
-    deformed_bond, deformed_bond_length = Geometry.bond_geometry(
-        nodes,
-        nlist,
-        deformed_coor,
-        deformed_bond,
-        deformed_bond_length,
-    )
+    Geometry.bond_geometry(nodes, nlist, deformed_coor, deformed_bond, deformed_bond_length)
     return datamanager
 end
 

diff --git a/src/Models/Pre_calculation/deformation_gradient.jl b/src/Models/Pre_calculation/deformation_gradient.jl
@@ -5,7 +5,7 @@
 module Deformation_Gradient
 using DataStructures: OrderedDict
 include("../../Support/geometry.jl")
-using .Geometry: compute_deformation_gradient
+using .Geometry: compute_deformation_gradients!
 export pre_calculation_name
 export init_model
 export compute
@@ -80,17 +80,18 @@ function compute(
     deformed_bond = datamanager.get_field("Deformed Bond Geometry", "NP1")
     deformation_gradient = datamanager.get_field("Deformation Gradient")
     inverse_shape_tensor = datamanager.get_field("Inverse Shape Tensor")
-
-    deformation_gradient = compute_deformation_gradient(
+    dof = datamanager.get_dof()
+    compute_deformation_gradients!(
+        deformation_gradient,
         nodes,
+        dof,
         nlist,
         volume,
         omega,
         bond_damage,
         deformed_bond,
         undeformed_bond,
         inverse_shape_tensor,
-        deformation_gradient,
     )
 
     return datamanager

diff --git a/src/Models/Pre_calculation/pre_bond_associated_correspondence.jl b/src/Models/Pre_calculation/pre_bond_associated_correspondence.jl
@@ -5,9 +5,10 @@
 module Pre_Bond_Associated_Correspondence
 using DataStructures: OrderedDict
 include("../../Support/geometry.jl")
-using .Geometry: calculate_bond_length, compute_weighted_deformation_gradient
+using .Geometry: compute_weighted_deformation_gradient
 using LoopVectorization
 using StaticArrays: @MVector
+
 export pre_calculation_name
 export init_model
 export compute
@@ -77,7 +78,9 @@ function init_model(
 
     accuracy_order = datamanager.get_accuracy_order()
     dim = qdim(accuracy_order, dof)
-
+    if "Q Vector" in datamanager.get_all_field_keys()
+        return datamanager
+    end
     datamanager.create_constant_free_size_field("Q Vector", Float64, (dim,), 1)
     datamanager.create_constant_free_size_field("Minv Matrix", Float64, (dim, dim))
     datamanager.create_constant_free_size_field("M temporary Matrix", Float64, (dim, dim))
@@ -124,8 +127,10 @@ function compute(
     Minv = datamanager.get_field("Minv Matrix")
     M = datamanager.get_field("M temporary Matrix")
 
-    weighted_volume =
-        compute_weighted_volume(nodes, nlist, volume, bond_damage, omega, weighted_volume)
+
+    compute_weighted_volume!(weighted_volume, nodes, nlist, volume, bond_damage, omega)
+
+
     compute_Lagrangian_gradient_weights(
         nodes,
         dof,
@@ -143,7 +148,6 @@ function compute(
     )
     compute_weighted_deformation_gradient(
         nodes,
-        dof,
         nlist,
         volume,
         gradient_weights,
@@ -156,22 +160,21 @@ end
 
 
 
-function compute_weighted_volume(
+function compute_weighted_volume!(
+    weighted_volume::Vector{Float64},
     nodes::Union{SubArray,Vector{Int64}},
     nlist::Union{Vector{Vector{Int64}},SubArray},
     volume::Vector{Float64},
     bond_damage::Vector{Vector{Float64}},
     omega::Vector{Vector{Float64}},
-    weighted_volume::Vector{Float64},
 )
     @views @inbounds @fastmath for iID in nodes
         weighted_volume[iID] = 0
         @views @inbounds @fastmath for jID in eachindex(nlist[iID])
-            @views weighted_volume[iID] +=
+            weighted_volume[iID] +=
                 bond_damage[iID][jID] * omega[iID][jID] * volume[nlist[iID][jID]]
         end
     end
-    return weighted_volume
 end
 
 
@@ -221,7 +224,6 @@ function calculate_Q(
 
 
     counter = 0
-    Q .= 1
     p = @MVector zeros(Int64, dof)
     for this_order = 1:accuracy_order
         for p[1] = this_order:-1:0
@@ -243,6 +245,7 @@ function calculate_Q(
 end
 
 function prod_Q(bond_geometry, horizon, p, Q)
+    Q = 1
     @inbounds @fastmath for m ∈ axes(p, 1)
         @views @inbounds @fastmath for pc = 1:p[m]
             Q *= bond_geometry[m] / horizon
@@ -270,6 +273,7 @@ function compute_Lagrangian_gradient_weights(
 
     for iID in nodes
         M .= 0
+        # bg = @view bond_geometry[iID]
         for (jID, nID) in enumerate(nlist[iID])
             Q = calculate_Q(accuracy_order, dof, bond_geometry[iID][jID], horizon[iID], Q)
             QTQ!(M, omega[iID][jID], bond_damage[iID][jID], volume[nID], Q)

diff --git a/src/Models/Pre_calculation/shape_tensor.jl b/src/Models/Pre_calculation/shape_tensor.jl
@@ -7,7 +7,7 @@ using DataStructures: OrderedDict
 include("../../Support/helpers.jl")
 using .Helpers: find_active_nodes
 include("../../Support/geometry.jl")
-using .Geometry
+using .Geometry: compute_shape_tensors
 export pre_calculation_name
 export init_model
 export compute
@@ -88,7 +88,7 @@ function compute(
     active_nodes = datamanager.get_field("Active Nodes")
     active_nodes = find_active_nodes(update_list, active_nodes, nodes)
 
-    shape_tensor, inverse_shape_tensor = Geometry.shape_tensor(
+    compute_shape_tensors(
         active_nodes,
         nlist,
         volume,