Skip to content

Commit 16b07ac

Browse files
CWillbergCWillberg
committed
code optimization
1 parent 558ba81 commit 16b07ac

File tree

8 files changed

+53
-45
lines changed

8 files changed

+53
-45
lines changed

‎src/Core/Solver/Matrix_Verlet.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -56,7 +56,7 @@ function init_solver(solver_options::Dict{Any,Any},
5656
thermal = "Thermal" in solver_options["Models"]
5757

5858
solver_options["dt"] = compute_crititical_time_step(block_nodes, mechanical, thermal)
59-
59+
#solver_options["dt"] = 1e-8
6060
solver_options["Number of Steps"] = Int(ceil((solver_options["Final Time"] -
6161
solver_options["Initial Time"]) /
6262
solver_options["dt"]))

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

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ using ........Helpers:
1212
find_local_neighbors, invert, rotate, determinant, smat,
1313
matrix_diff!
1414
using ........Geometry:
15-
compute_strain,
15+
compute_strain!,
1616
compute_bond_level_rotation_tensor,
1717
compute_bond_level_deformation_gradient
1818
using ....Pre_Calculation.Pre_Bond_Associated_Correspondence: compute_weighted_volume!
@@ -269,8 +269,8 @@ function compute_bond_strain(nodes,
269269
strain_N,
270270
strain_increment)
271271
for iID in nodes
272-
compute_strain(eachindex(nlist[iID]), deformation_gradient[iID],
273-
strain_NP1[iID])
272+
compute_strain!(eachindex(nlist[iID]), deformation_gradient[iID],
273+
strain_NP1[iID])
274274
matrix_diff!(strain_increment[iID],
275275
eachindex(nlist[iID]),
276276
strain_NP1[iID],

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

Lines changed: 9 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ include("./Bond_Associated_Correspondence.jl")
2222
using .Bond_Associated_Correspondence
2323
using ....Material_Basis: compute_Piola_Kirchhoff_stress!
2424
using .......Helpers: invert, rotate, determinant, smat, matrix_diff!, fast_mul!, mat_mul!
25-
using .......Geometry: compute_strain
25+
using .......Geometry: compute_strain!
2626

2727
export init_model
2828
export material_name
@@ -161,7 +161,7 @@ function compute_correspondence_model(nodes::AbstractVector{Int64},
161161
time::Float64,
162162
dt::Float64)
163163
rotation::Bool = Data_Manager.get_rotation()
164-
dof = Data_Manager.get_dof()
164+
dof::Int64 = Data_Manager.get_dof()
165165
deformation_gradient::NodeTensorField{Float64} = Data_Manager.get_field("Deformation Gradient")
166166
bond_force::BondVectorState{Float64} = Data_Manager.get_field("Bond Forces")
167167
bond_damage::BondScalarState{Float64} = Data_Manager.get_bond_damage("NP1")
@@ -173,8 +173,10 @@ function compute_correspondence_model(nodes::AbstractVector{Int64},
173173
stress_N::NodeTensorField{Float64} = Data_Manager.get_field("Cauchy Stress", "N")
174174
stress_NP1::NodeTensorField{Float64} = Data_Manager.get_field("Cauchy Stress", "NP1")
175175
strain_increment::NodeTensorField{Float64} = Data_Manager.get_field("Strain Increment")
176-
compute_strain(nodes, deformation_gradient, strain_NP1)
177-
matrix_diff!(strain_increment, nodes, strain_NP1, strain_N)
176+
177+
@timeit "compute strain" compute_strain!(nodes, deformation_gradient, strain_NP1)
178+
@timeit "compute matrix diff" matrix_diff!(strain_increment, nodes, strain_NP1,
179+
strain_N)
178180

179181
if rotation
180182
rotation_tensor::NodeTensorField{Float64} = Data_Manager.get_field("Rotation Tensor")
@@ -184,7 +186,7 @@ function compute_correspondence_model(nodes::AbstractVector{Int64},
184186

185187
# material_models = split(material_parameter["Material Model"], "+")
186188
# material_models = map(r -> strip(r), material_models)
187-
@timeit "Calculate material" begin
189+
@timeit "compute material" begin
188190
for material_model in Data_Manager.get_analysis_model("Correspondence Model",
189191
block)
190192
mod::Module = Data_Manager.get_model_module(material_model)
@@ -201,9 +203,9 @@ function compute_correspondence_model(nodes::AbstractVector{Int64},
201203
end
202204

203205
if rotation
204-
stress_NP1 = rotate(nodes, stress_NP1, rotation_tensor, true)
206+
rotate(nodes, stress_NP1, rotation_tensor, true)
205207
end
206-
@timeit "Compute bond force" calculate_bond_force!(nodes,
208+
@timeit "compute bond force" calculate_bond_force!(nodes,
207209
dof,
208210
deformation_gradient,
209211
undeformed_bond,

‎src/Models/Material/Material_Models/Correspondence/Correspondence_Elastic.jl‎

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,6 @@
33
# SPDX-License-Identifier: BSD-3-Clause
44

55
module Correspondence_Elastic
6-
76
using .......Data_Manager
87
using .....Material_Basis: get_Hooke_matrix
98
using .......Helpers: get_fourth_order, fast_mul!, get_mapping
@@ -110,7 +109,14 @@ function compute_stresses(nodes::AbstractVector{Int64},
110109
strain_increment::NodeTensorField{Float64},
111110
stress_N::NodeTensorField{Float64},
112111
stress_NP1::NodeTensorField{Float64})
113-
mapping = get_mapping(dof)
112+
mapping = if dof == 2
113+
get_mapping(2)::SMatrix{3,2,Int64,6}
114+
elseif dof == 3
115+
get_mapping(3)::SMatrix{6,2,Int64,12}
116+
else
117+
get_mapping(dof)
118+
end
119+
114120
hooke_matrix::NodeTensorField{Float64} = Data_Manager.get_field("Elasticity Matrix")
115121

116122
n_map = size(mapping, 1)

‎src/Models/Material/Material_Models/Correspondence/Correspondence_matrix_based.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -632,7 +632,7 @@ function init_matrix()
632632
inverse_shape_tensor,
633633
nlist,
634634
volume,
635-
bond_geometry_N,
635+
bond_geometry,
636636
bond_damage,
637637
omega)
638638
end

‎src/Support/Geometry.jl‎

Lines changed: 23 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,7 @@ using StaticArrays
99
using Rotations
1010
using LoopVectorization: @avx, @fastmath, @simd
1111
using ...Data_Manager
12-
using ...Helpers: invert, smat, mat_mul_transpose_mat!
12+
using ...Helpers: invert, smat
1313

1414
export bond_geometry!
1515
export compute_shape_tensors!
@@ -340,7 +340,9 @@ function deformation_gradient_decomposition(nodes::Union{Base.OneTo{Int64},Vecto
340340
end
341341

342342
"""
343-
function compute_strain(nodes::Union{Base.OneTo{Int64},Vector{Int64}, SubArray}, deformation_gradient, strain)
343+
function compute_strain!(nodes::AbstractVector{Int64},
344+
deformation_gradient::NodeTensorField{Float64},
345+
strain::NodeTensorField{Float64})
344346
345347
Calculate strains for specified nodes based on deformation gradients.
346348
@@ -355,17 +357,26 @@ Calculate strains for specified nodes based on deformation gradients.
355357
This function iterates over the specified nodes and computes strain at each node using the given deformation gradients.
356358
357359
"""
358-
function compute_strain(nodes::AbstractVector{Int64},
359-
deformation_gradient::NodeTensorField{Float64},
360-
strain::NodeTensorField{Float64})
361-
# https://en.wikipedia.org/wiki/Strain_(mechanics)
360+
function compute_strain!(nodes::AbstractVector{Int64},
361+
deformation_gradient::NodeTensorField{Float64},
362+
strain::NodeTensorField{Float64})
363+
# Green-Lagrange strain: E = 0.5*(F^T*F - I)
364+
362365
for iID in nodes
363-
@views mat_mul_transpose_mat!(strain[iID, :, :],
364-
deformation_gradient[iID, :, :],
365-
deformation_gradient[iID, :, :])
366-
if !all(deformation_gradient[iID, :, :] .== 0.0)
367-
@inbounds @fastmath for m in axes(strain, 2)
368-
strain[iID, m, m] -= 0.5
366+
F = @view deformation_gradient[iID, :, :]
367+
E = @view strain[iID, :, :]
368+
369+
@inbounds @fastmath for m in axes(F, 1)
370+
for n in m:size(F, 2) # n >= m
371+
Cmn = m == n ? -1.0 : 0.0
372+
@inbounds @fastmath for k in axes(F, 2)
373+
Cmn += F[k, m] * F[k, n]
374+
end
375+
E[m, n] = 0.5 * Cmn
376+
# Symmetry: E[n, m] = E[m, n]
377+
if m != n
378+
E[n, m] = E[m, n]
379+
end
369380
end
370381
end
371382
end

‎src/Support/Helpers.jl‎

Lines changed: 4 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@ using LoopVectorization
1515
using Unitful: ustrip
1616
using CDDLib
1717
using Polyhedra: MixedMatHRep, polyhedron, vrep, hrep
18-
18+
using Base.Threads: @threads
1919
using ..Data_Manager
2020

2121
export qdim
@@ -43,7 +43,6 @@ export fastdot
4343
export fast_mul!
4444
export mat_mul!
4545
export get_mapping
46-
export mat_mul_transpose_mat!
4746
export get_ring
4847
export get_hexagon
4948
export nearest_point_id
@@ -371,9 +370,8 @@ end
371370
function matrix_diff!(s3::AbstractArray{Float64,3}, nodes::AbstractVector{Int64},
372371
s2::AbstractArray{Float64,3},
373372
s1::AbstractArray{Float64,3})
374-
@views @inbounds for iID in nodes
375-
@inbounds @fastmath @views for m in axes(s1, 2),
376-
n in axes(s1, 3)
373+
@inbounds for iID in nodes
374+
@inbounds @fastmath for m in axes(s1, 2), n in axes(s1, 3)
377375
s3[iID, m, n] = s2[iID, m, n] - s1[iID, m, n]
378376
end
379377
end
@@ -413,16 +411,7 @@ function mat_mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T},
413411
C[m, n] = Cmn
414412
end
415413
end
416-
function mat_mul_transpose_mat!(C::AbstractMatrix{T}, A::AbstractMatrix{T},
417-
B::AbstractMatrix{T}) where {T<:Float64}
418-
@inbounds @fastmath for m in axes(A, 1), n in axes(B, 2)
419-
Cmn = zero(T)
420-
@inbounds @fastmath for k in axes(A, 2)
421-
Cmn += A[k, m] * B[k, n]
422-
end
423-
C[m, n] = 0.5 * Cmn
424-
end
425-
end
414+
426415
function add_in_place!(C::AbstractMatrix{T},
427416
A::Vector{Vector{T}},
428417
B::Vector{Vector{T}},

‎test/unit_tests/Support/ut_geometry.jl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -451,7 +451,7 @@ end
451451
undeformed_bond,
452452
undeformed_bond,
453453
inverse_shape_tensor)
454-
PeriLab.Geometry.compute_strain(nodes, deformation_gradient, strain)
454+
PeriLab.Geometry.compute_strain!(nodes, deformation_gradient, strain)
455455

456456
for i in 1:nnodes
457457
@test strain[i, 1, 1] == 0
@@ -473,9 +473,9 @@ end
473473
deformation_gradient_3D[1, 3, 3] = 3.0
474474
strain_3D = test_data_manager.create_constant_node_tensor_field("Strain_3D", Float64,
475475
3)
476-
PeriLab.Geometry.compute_strain(view(nodes, eachindex(nodes)),
477-
deformation_gradient_3D,
478-
strain_3D)
476+
PeriLab.Geometry.compute_strain!(view(nodes, eachindex(nodes)),
477+
deformation_gradient_3D,
478+
strain_3D)
479479
identity = zeros(3, 3)
480480
identity[1, 1] = 1
481481
identity[2, 2] = 1

0 commit comments

Comments
 (0)