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

Project.toml

@@ -34,7 +34,6 @@ Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
-Tensorial = "98f94333-fa9f-48a9-ad80-1c66397b2b38"
 Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
@@ -72,7 +71,6 @@ Rotations = "^1"
 StaticArrays = "^1"
 Statistics = "^1"
 TensorOperations = "^4, 5"
-Tensorial = "0.16, 0.17, 0.18"
 Tensors = "^1"
 Test = "^1"
 TestSetExtensions = "^3"

src/Compute/compute_field_values.jl

@@ -2,7 +2,7 @@
 #
 # SPDX-License-Identifier: BSD-3-Clause
 
-using .Helpers: find_active_nodes, get_active_update_nodes
+using .Helpers: find_active_nodes, get_active_update_nodes, add_in_place!
 using StaticArrays: MMatrix, SMatrix
 using .Helpers: invert
 include("../Models/Material/material_basis.jl")
@@ -47,7 +47,8 @@ function get_partial_stresses(datamanager::Module, nodes::Union{SubArray,Vector{
     stress = datamanager.get_field("Cauchy Stress", "NP1")
 
     for iID in nodes
-        stress[iID, :, :] .+= bond_forces[iID]' * undeformed_bond[iID] .* volume[iID]
+        str = @view stress[iID, :, :]
+        add_in_place!(str, bond_forces[iID], undeformed_bond[iID], volume[iID])
     end
     return datamanager
 end

src/Core/data_manager.jl

@@ -41,7 +41,6 @@ export get_num_responder
 export get_max_rank
 export get_cancel
 export get_damage_models
-export get_material_model
 export get_output_frequency
 export get_rotation
 export get_element_rotation
@@ -56,7 +55,6 @@ export set_inverse_nlist
 export set_fem
 export set_glob_to_loc
 export set_damage_models
-export set_material_model
 export set_model_module
 export set_num_controller
 export set_nset
@@ -453,6 +451,7 @@ function create_field(
                 fields[vartype][name] = fill(value, (data["nnodes"], dof, dof))
             else
                 fields[vartype][name] = fill(value, data["nnodes"], dof)
+                # fields[vartype][name] = [fill(value,dof) for j=1:data["nnodes"]]
             end
         end
     elseif bond_or_node == "Bond_Field"
@@ -463,7 +462,8 @@ function create_field(
             if VectorOrArray == "Matrix"
                 fields[vartype][name] = [fill(value, (n, dof, dof)) for n in nBonds]
             else
-                fields[vartype][name] = [fill(value, (n, dof)) for n in nBonds]
+                # fields[vartype][name] = [fill(value, (n, dof)) for n in nBonds]
+                fields[vartype][name] = [[fill(value, dof) for j = 1:n] for n in nBonds]
             end
         end
     elseif bond_or_node == "Element_Field"
@@ -990,18 +990,6 @@ function get_damage_models()
     return data["damage_models"]
 end
 
-"""
-    get_material_models()
-
-This function returns the `material_models` variable.
-
-# Returns
-- `material_models`::Any: The value of the `material_models` variable.
-"""
-function get_material_models()
-    return data["material_models"]
-end
-
 """
     loc_to_glob(range::UnitRange{Int64})
 
@@ -1583,8 +1571,13 @@ function switch_NP1_to_N(nstep::Int64)
 end
 
 function switch_nodes!(field_N, field_NP1, NP1)
-    copyto!(field_N, field_NP1)
-    fill!(field_NP1, data["field_types"][NP1](0))
+    if field_NP1[1] isa AbstractVector
+        copyto!.(field_N, field_NP1)
+        fill!.(field_NP1, data["field_types"][NP1](0))
+    else
+        copyto!(field_N, field_NP1)
+        fill!(field_NP1, data["field_types"][NP1](0))
+    end
 end
 function switch_bonds!(field_N, field_NP1, NP1)
 

src/IO/mesh_data.jl

@@ -169,7 +169,7 @@ function create_and_distribute_bond_norm(
     datamanager::Module,
     nlist_filtered_ids::Vector{Vector{Int64}},
     distribution::Vector{Vector{Int64}},
-    bond_norm::Vector{Matrix{Float64}},
+    bond_norm::Vector{Any},
     dof::Int64,
 )
     bond_norm = send_value(comm, 0, bond_norm)
@@ -181,7 +181,7 @@ function create_and_distribute_bond_norm(
         distribution,
         true,
     )
-    bond_norm_field .= bond_norm
+    copyto!.(bond_norm_field, bond_norm)
 end
 
 """
@@ -1705,11 +1705,9 @@ function apply_bond_filters(
         end
         if contact_enabled
             nlist_filtered_ids = fill(Vector{Int64}([]), nnodes)
-            bond_norm = Matrix{Float64}[]
+            bond_norm = []
             for iID = 1:nnodes
-                # bond_norm[iID] = fill(fill(1, dof), length(nlist[iID]))
-                append!(bond_norm, [Matrix{Float64}(undef, length(nlist[iID]), dof)])
-                fill!(bond_norm[end], Float64(1))
+                push!(bond_norm, [fill(1.0, dof) for n = 1:length(nlist[iID])])
             end
         end
 
@@ -1728,7 +1726,7 @@ function apply_bond_filters(
                     )
                     nlist_filtered_ids[iID] = indices
                     for jID in indices
-                        bond_norm[iID][jID, :] .= normal
+                        bond_norm[iID][jID] .= normal
                     end
                 else
                     nlist[iID] = nlist[iID][filter_flag[iID]]

src/MPI_communication/MPI_communication.jl

@@ -234,7 +234,7 @@ function synch_controller_bonds_to_responder(comm::MPI.Comm, overlapnodes, array
                     @views send_msg = array[iID]
                 else
                     #TODO: Check if we can remove the [:,:]
-                    @views send_msg = array[iID][:, :]
+                    @views send_msg = mapreduce(permutedims, vcat, array[iID])
                 end
                 MPI.Isend(send_msg, comm; dest = jcore - 1, tag = 0)
                 # @debug "Sending   $rank -> $(jcore-1)"
@@ -245,15 +245,15 @@ function synch_controller_bonds_to_responder(comm::MPI.Comm, overlapnodes, array
                 if dof == 1
                     @views recv_msg = similar(array[iID])
                 else
-                    @views recv_msg = similar(array[iID][:, :])
+                    @views recv_msg = similar(mapreduce(permutedims, vcat, array[iID]))
                 end
                 MPI.Recv!(recv_msg, comm; source = jcore - 1, tag = 0)
                 # @debug "Receiving $(jcore-1) -> $rank"
                 recv_msg = reshape(recv_msg, :, 2)
                 if dof == 1
                     array[iID] = recv_msg
                 else
-                    array[iID][:, :] = recv_msg
+                    array[iID] = Vector{eltype(recv_msg)}[eachcol(recv_msg)...]
                 end
             end
         end
@@ -274,11 +274,11 @@ Split a vector into a vector of matrices
 - `result::Vector`: The result vector
 """
 function split_vector(input, row_nums, dof)
-    result = Vector{Matrix{eltype(input)}}()
+    result = Vector{Vector{Vector{eltype(input)}}}()
     start = firstindex(input)
-    for len in row_nums
-        push!(result, reshape(input[start:(start+len-1)], (Int64(len / dof), dof)))
-        start += len
+    for (i, len) in enumerate(row_nums)
+        push!(result, [input[start+(j-1)*dof:start-1+j*dof] for j = 1:len])
+        start += dof * len
     end
     result
 end
@@ -317,7 +317,7 @@ function synch_controller_bonds_to_responder_flattened(
             @views send_indices = overlapnodes[rank+1][jcore]["Controller"]
             # @debug "Sending $rank -> $(jcore-1)"
             MPI.Isend(
-                vcat([reshape(array[i], :, 1) for i in send_indices]...),
+                vcat(vcat(array...)[send_indices]...),
                 comm;
                 dest = jcore - 1,
                 tag = 0,
@@ -326,7 +326,7 @@ function synch_controller_bonds_to_responder_flattened(
         if !isempty(overlapnodes[rank+1][jcore]["Responder"])
             @views recv_indices = overlapnodes[rank+1][jcore]["Responder"]
             row_nums = [length(subarr) for subarr in array[recv_indices]]
-            @views recv_msg = zeros(sum(row_nums))
+            @views recv_msg = zeros(sum(row_nums * dof))
             # @debug "Receiving $(jcore-1) -> $rank"
             MPI.Recv!(recv_msg, comm; source = jcore - 1, tag = 0)
             recv_msg = split_vector(recv_msg, row_nums, dof)

src/Models/Damage/Energy_release.jl

@@ -8,7 +8,7 @@ include("../../Support/geometry.jl")
 include("../../Support/helpers.jl")
 using .Material
 using .Geometry
-using .Helpers: rotate, fastdot, fastsubtract!
+using .Helpers: rotate, fastdot, sub_in_place!
 using LinearAlgebra
 using StaticArrays
 export compute_model
@@ -77,7 +77,7 @@ function compute_model(
     critical_energy =
         critical_field ? datamanager.get_field("Critical_Value") :
         damage_parameter["Critical Value"]
-    quad_horizon = datamanager.get_field("Quad Horizon")
+    quad_horizons = datamanager.get_field("Quad Horizon")
     inverse_nlist = datamanager.get_inverse_nlist()
     dependend_value::Bool = false
     if haskey(damage_parameter, "Temperature dependend")
@@ -109,17 +109,26 @@ function compute_model(
 
     bond_energy::Float64 = 0.0
     norm_displacement::Float64 = 0.0
+    product::Float64 = 0.0
     x_vector::Vector{Float64} = @SVector zeros(Float64, dof)
     x_vector[1] = 1.0
 
-    neighbor_bond_force::Vector{Float64} = @SVector zeros(Float64, dof)
-    projected_force::Vector{Float64} = @SVector zeros(Float64, dof)
+    bond_force::Vector{Float64} = zeros(Float64, dof)
+    neighbor_bond_force::Vector{Float64} = zeros(Float64, dof)
+    bond_force_delta::Vector{Float64} = zeros(Float64, dof)
+    projected_force::Vector{Float64} = zeros(Float64, dof)
+    relative_displacement::Vector{Float64} = zeros(Float64, dof)
+
+    sub_in_place!(bond_displacements, deformed_bond, undeformed_bond)
 
-    fastsubtract!(bond_displacements, deformed_bond, undeformed_bond)
     for iID in nodes
-        @views nlist_temp = nlist[iID]
+        nlist_temp = nlist[iID]
+        bond_displacement = bond_displacements[iID]
+        bond_force_vec = bond_forces[iID]
+        quad_horizon = quad_horizons[iID]
         for jID in eachindex(nlist_temp)
-            @views relative_displacement = bond_displacements[iID][jID, :]
+
+            relative_displacement .= bond_displacement[jID]
             norm_displacement = fastdot(relative_displacement, relative_displacement)
             if norm_displacement == 0 || (
                 tension &&
@@ -132,20 +141,20 @@ function compute_model(
 
             # check if the bond also exist at other node, due to different horizons
             try
-                @views neighbor_bond_force .=
-                    bond_forces[neighborID][inverse_nlist[neighborID][iID], :]
+                neighbor_bond_force .=
+                    bond_forces[neighborID][inverse_nlist[neighborID][iID]]
             catch e
                 # Handle the case when the key doesn't exist
             end
 
-            @views projected_force .=
-                fastdot(
-                    (bond_forces[iID][jID, :] - neighbor_bond_force),
-                    relative_displacement,
-                ) / (norm_displacement) .* relative_displacement
+            bond_force .= bond_force_vec[jID]
+            bond_force_delta .= bond_force .- neighbor_bond_force
+
 
-            @views bond_energy =
-                0.25 * fastdot((abs.(projected_force)), (abs.(relative_displacement)))
+            product = fastdot(bond_force_delta, relative_displacement)
+            mul!(projected_force, product / norm_displacement, relative_displacement)
+            product = fastdot(projected_force, relative_displacement, true)
+            bond_energy = 0.25 * product
 
             if dependend_value
                 critical_energy =
@@ -160,10 +169,10 @@ function compute_model(
 
             if aniso_damage
                 if all(rotation_tensor[iID, :, :] .== 0)
-                    @views rotated_bond = deformed_bond[iID][jID, :]
+                    @views rotated_bond = deformed_bond[iID][jID]
                 else
                     @views rotated_bond =
-                        rotation_tensor[iID, :, :]' * deformed_bond[iID][jID, :]
+                        rotation_tensor[iID, :, :]' * deformed_bond[iID][jID]
                 end
                 # Compute bond_norm for all components at once
                 @views bond_norm_all = abs.(rotated_bond) ./ deformed_bond_length[iID][jID]
@@ -207,7 +216,8 @@ function compute_model(
                 # end
 
             else
-                if bond_energy > crit_energy * quad_horizon[iID]
+                product = crit_energy * quad_horizon
+                if bond_energy > product
                     bond_damage[iID][jID] = 0.0
                     update_list[iID] = true
                 end

src/Models/Material/Material_Factory.jl

@@ -43,7 +43,7 @@ function init_fields(datamanager::Module)
     datamanager.create_constant_node_field("Acceleration", Float64, dof)
     datamanager.create_node_field("Velocity", Float64, dof)
     datamanager.create_constant_bond_field("Bond Forces", Float64, dof)
-    dof = datamanager.get_dof()
+    datamanager.create_constant_bond_field("Temporary Bond Field", Float64, 1)
     deformed_coorN, deformed_coorNP1 =
         datamanager.create_node_field("Deformed Coordinates", Float64, dof)
     deformed_coorN = copy(datamanager.get_field("Coordinates"))
@@ -111,11 +111,8 @@ function init_model(datamanager::Module, nodes::Union{SubArray,Vector{Int64}}, b
         for iID in nodes
             if length(nlist_filtered_ids[iID]) != 0
                 for neighborID in nlist_filtered_ids[iID]
-                    bond_norm[iID][neighborID, :] .*= sign(
-                        dot(
-                            (bond_geometry[iID][neighborID, :]),
-                            bond_norm[iID][neighborID, :],
-                        ),
+                    bond_norm[iID][neighborID] .*= sign(
+                        dot((bond_geometry[iID][neighborID]), bond_norm[iID][neighborID]),
                     )
                 end
             end

src/Models/Material/Material_Models/BondBased/Bondbased_Elastic.jl

@@ -147,14 +147,14 @@ function compute_bb_force!(
     deformed_bond,
 )
     @inbounds @fastmath for i ∈ axes(bond_force, 1)
-        @inbounds @fastmath for j ∈ axes(bond_force, 2)
-            bond_force[i, j] =
+        @inbounds @fastmath for j ∈ axes(bond_force[i], 1)
+            bond_force[i][j] =
                 (
                     constant *
                     bond_damage[i] *
                     (deformed_bond_length[i] - undeformed_bond_length[i]) /
                     undeformed_bond_length[i]
-                ) * deformed_bond[i, j] / deformed_bond_length[i]
+                ) * deformed_bond[i][j] / deformed_bond_length[i]
         end
     end
 end

src/Models/Material/Material_Models/Correspondence/Bond_Associated_Correspondence.jl

@@ -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
@@ -280,7 +279,7 @@ function compute_stress_integral(
     bond_damage::Vector{Vector{Float64}},
     volume::Vector{Float64},
     weighted_volume::Vector{Float64},
-    bond_geometry::Vector{Matrix{Float64}},
+    bond_geometry::Vector{Vector{Vector{Float64}}},
     bond_length::Vector{Vector{Float64}},
     bond_stresses::Vector{Array{Float64,3}},
     deformation_gradient::Vector{Array{Float64,3}},
@@ -300,7 +299,7 @@ function compute_stress_integral(
                 continue
             end
             @views temp =
-                (one - bond_geometry[iID][jID, :] * bond_geometry[iID][jID, :]') ./
+                (one - bond_geometry[iID][jID] * bond_geometry[iID][jID]') ./
                 (bond_length[iID][jID] * bond_length[iID][jID])
 
             @views factor =
@@ -430,13 +429,13 @@ function compute_bond_forces(
 
             # bond_forces[iID][jID, :] =
             #     integral_nodal_stress[iID, :, :] * gradient_weights[iID][jID, :]
-            bf = @view bond_forces[iID][jID, :]
-            gw = @view 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, :] +=
+            @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, :]
+                bond_stress[iID][jID, :, :] * bond_geometry[iID][jID]
 
 
             #