From 4046f5288d7f01c36df164aa82bf41e3c43fe1b2 Mon Sep 17 00:00:00 2001 From: will_cr Date: Mon, 28 Oct 2024 14:29:40 +0100 Subject: [PATCH 1/6] memory optimization --- .../Bond_Associated_Correspondence.jl | 11 +++-------- src/Support/geometry.jl | 15 ++++++++++----- src/Support/helpers.jl | 11 ++++++++++- .../ut_Bond_Associated_Correspondence.jl | 8 ++++---- 4 files changed, 27 insertions(+), 18 deletions(-) diff --git a/src/Models/Material/Material_Models/Correspondence/Bond_Associated_Correspondence.jl b/src/Models/Material/Material_Models/Correspondence/Bond_Associated_Correspondence.jl index 51a07b11..60ef363d 100644 --- a/src/Models/Material/Material_Models/Correspondence/Bond_Associated_Correspondence.jl +++ b/src/Models/Material/Material_Models/Correspondence/Bond_Associated_Correspondence.jl @@ -175,7 +175,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,12 +336,8 @@ 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], @@ -349,7 +345,6 @@ function compute_bond_strain( ) end - return strain_NP1, strain_increment end diff --git a/src/Support/geometry.jl b/src/Support/geometry.jl index 92babb29..8d517996 100644 --- a/src/Support/geometry.jl +++ b/src/Support/geometry.jl @@ -8,7 +8,7 @@ using TensorOperations using Combinatorics: levicivita using Rotations include("helpers.jl") -using .Helpers: invert, smat +using .Helpers: invert, smat, mat_mul_transpose_mat! export bond_geometry export shape_tensor @@ -345,13 +345,18 @@ function compute_strain( deformation_gradient::Union{Array{Float64,3},SubArray}, strain::Union{Array{Float64,3},SubArray}, ) - # https://en.wikipedia.org/wiki/Strain_(mechanics) for iID in nodes - def_grad = smat(deformation_gradient[iID, :, :]) - strain[iID, :, :] = 0.5 .* (def_grad' * def_grad - I) + @views mat_mul_transpose_mat!( + strain[iID, :, :], + deformation_gradient[iID, :, :], + deformation_gradient[iID, :, :], + ) + @inbounds @fastmath for m ∈ axes(strain, 2) + strain[iID, m, m] -= 0.5 + end end - return strain + end """ diff --git a/src/Support/helpers.jl b/src/Support/helpers.jl index 929032e3..f9ada13d 100644 --- a/src/Support/helpers.jl +++ b/src/Support/helpers.jl @@ -35,6 +35,7 @@ export fastdot export fast_mul! export mat_mul! export get_mapping +export mat_mul_transpose_mat! @@ -78,7 +79,15 @@ function mat_mul!(C, A, B) C[m, n] = Cmn end end - +function mat_mul_transpose_mat!(C, A, B) + @inbounds @fastmath for m ∈ axes(A, 1), n ∈ axes(B, 2) + Cmn = zero(eltype(C)) + @inbounds @fastmath for k ∈ axes(A, 2) + Cmn += A[k, m] * B[k, n] + end + C[m, n] = 0.5 * Cmn + end +end function get_MMatrix(len::Int64) diff --git a/test/unit_tests/Models/Material/Material_Models/Correspondence/ut_Bond_Associated_Correspondence.jl b/test/unit_tests/Models/Material/Material_Models/Correspondence/ut_Bond_Associated_Correspondence.jl index 0ee7ca7a..edde1e94 100644 --- a/test/unit_tests/Models/Material/Material_Models/Correspondence/ut_Bond_Associated_Correspondence.jl +++ b/test/unit_tests/Models/Material/Material_Models/Correspondence/ut_Bond_Associated_Correspondence.jl @@ -6,8 +6,8 @@ include( "../../../../../../src/Models/Material/Material_Models/Correspondence/Bond_Associated_Correspondence.jl", ) using Test -#include("../../../../../../src/PeriLab.jl") -#using .PeriLab +include("../../../../../../src/PeriLab.jl") +using .PeriLab @testset "ut_ba_correspondence_name" begin @@ -49,7 +49,7 @@ end strainN[1][1, :, :] = [-1 -1; -1 -1] strainN[2][1, :, :] = [0.5 0; 0 0.5] - strain, strain_inc = Bond_Associated_Correspondence.compute_bond_strain( + Bond_Associated_Correspondence.compute_bond_strain( nodes, nlist, deformation_gradient, @@ -61,7 +61,7 @@ end @test strain[1][1, :, :] == [0 0; 0 0] @test strain[2][1, :, :] == [-0.5 0; 0 0.5] @test strain_inc[1][1, :, :] == [1 1; 1 1] - @test strain_inc[2][1, :, :] == [-1.0 0; 0 0.0] + @test strain_inc[2][1, :, :] == [-1.0 0.0; 0.0 0.0] end @testset "ut_update_Green_Langrange_strain" begin From e38b9d1b1303d4a1f47657ca65b6cd9b22760115 Mon Sep 17 00:00:00 2001 From: will_cr Date: Tue, 29 Oct 2024 17:29:38 +0100 Subject: [PATCH 2/6] memory optimization --- src/IO/mesh_data.jl | 4 +- .../Bond_Associated_Correspondence.jl | 10 +++-- .../Correspondence/Correspondence.jl | 3 +- .../Correspondence/Correspondence_Elastic.jl | 11 ++--- .../Pre_calculation/bond_deformation.jl | 8 +--- .../pre_bond_associated_correspondence.jl | 34 ++++++---------- src/Support/geometry.jl | 39 +++++++++++------- test/unit_tests/Core/Solver/ut_Verlet.jl | 2 +- test/unit_tests/IO/ut_mesh_data.jl | 2 +- .../ut_pre_bond_associated_correspondence.jl | 4 +- test/unit_tests/Support/ut_geometry.jl | 40 ++++++++++--------- 11 files changed, 77 insertions(+), 80 deletions(-) diff --git a/src/IO/mesh_data.jl b/src/IO/mesh_data.jl index 0878c916..d53d5643 100644 --- 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 index 60ef363d..3cbe6442 100644 --- a/src/Models/Material/Material_Models/Correspondence/Bond_Associated_Correspondence.jl +++ b/src/Models/Material/Material_Models/Correspondence/Bond_Associated_Correspondence.jl @@ -432,10 +432,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 index 83c83e62..c62805a0 100644 --- 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 index b6e713e5..e6b349f9 100644 --- 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 index 32e7e43b..fbce4091 100644 --- 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/pre_bond_associated_correspondence.jl b/src/Models/Pre_calculation/pre_bond_associated_correspondence.jl index f63f7b3a..e7472830 100644 --- a/src/Models/Pre_calculation/pre_bond_associated_correspondence.jl +++ b/src/Models/Pre_calculation/pre_bond_associated_correspondence.jl @@ -8,6 +8,7 @@ include("../../Support/geometry.jl") using .Geometry: calculate_bond_length, 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(nodes, nlist, volume, bond_damage, omega, weighted_volume) + + compute_Lagrangian_gradient_weights( nodes, dof, @@ -143,7 +148,6 @@ function compute( ) compute_weighted_deformation_gradient( nodes, - dof, nlist, volume, gradient_weights, @@ -171,7 +175,6 @@ function compute_weighted_volume( bond_damage[iID][jID] * omega[iID][jID] * volume[nlist[iID][jID]] end end - return weighted_volume end @@ -214,14 +217,13 @@ end function calculate_Q( accuracy_order::Int64, dof::Int64, - bond_geometry::Vector{Float64}, + bond_geometry, horizon::Union{Int64,Float64}, Q::Vector{Float64}, ) 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,14 +273,9 @@ 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, - ) + @views Q = calculate_Q(accuracy_order, dof, bg[jID, :], horizon[iID], Q) QTQ!(M, omega[iID][jID], bond_damage[iID][jID], volume[nID], Q) end @@ -287,13 +285,7 @@ function compute_Lagrangian_gradient_weights( ) for jID in eachindex(nlist[iID]) - Q = calculate_Q( - accuracy_order, - dof, - bond_geometry[iID][jID, :], - horizon[iID], - Q, - ) + @views Q = calculate_Q(accuracy_order, dof, bg[jID, :], horizon[iID], Q) # this comes from Eq(19) in 10.1007/s40571-019-00266-9 # or example 1 in https://arxiv.org/pdf/2004.11477 # Eq (3) flowing diff --git a/src/Support/geometry.jl b/src/Support/geometry.jl index 8d517996..97050d24 100644 --- a/src/Support/geometry.jl +++ b/src/Support/geometry.jl @@ -6,6 +6,7 @@ module Geometry using LinearAlgebra using TensorOperations using Combinatorics: levicivita +using StaticArrays using Rotations include("helpers.jl") using .Helpers: invert, smat, mat_mul_transpose_mat! @@ -43,6 +44,7 @@ Calculate bond geometries between nodes based on their coordinates. undeformed_bond(nodes, dof, nlist, coor, undeformed_bond) """ + function bond_geometry( nodes::Union{SubArray,Vector{Int64}}, nlist::Vector{Vector{Int64}}, @@ -52,29 +54,37 @@ function bond_geometry( ) for iID in nodes - undeformed_bond[iID], undeformed_bond_length[iID] = - calculate_bond_length(iID, coor, nlist[iID]) - end - if any(any.(x -> x == 0, [lengths for lengths in undeformed_bond_length[nodes]])) - @error "Identical point coordinates with no distance" - return nothing + calculate_bond_length( + iID, + coor, + nlist[iID], + undeformed_bond[iID], + undeformed_bond_length[iID], + ) end - return undeformed_bond, undeformed_bond_length + end + function calculate_bond_length( iID::Int64, coor::Union{SubArray,Matrix{Float64},Matrix{Int64}}, nlist::Vector{Int64}, + bond_vectors, + bond_norm, ) + @inbounds @fastmath @views for m ∈ axes(nlist, 1), cmn in zero(eltype(coor)) - @views bond_vectors = coor[nlist, :] .- coor[iID, :]' - # distances = sqrt.(sum(bond_vectors .^ 2, dims=2))[:] - - # Check for identical point coordinates - return bond_vectors, norm.(eachrow(bond_vectors)) + @inbounds @fastmath @views for n ∈ axes(coor, 2) + bond_vectors[m, n] = coor[nlist[m], n] - coor[iID, n] + cmn += bond_vectors[m, n] * bond_vectors[m, n] + end + if cmn == 0 + @error "Identical point coordinates with no distance" + end + bond_norm[m] = sqrt(cmn) + end end - """ shape_tensor(nodes::Union{SubArray, Vector{Int64}}, nlist, volume, omega, bond_damage, undeformed_bond, shape_tensor, inverse_shape_tensor) @@ -281,7 +291,6 @@ end function compute_weighted_deformation_gradient( nodes::Union{SubArray,Vector{Int64}}, - dof::Int64, nlist, volume, gradient_weight, @@ -407,7 +416,7 @@ function compute_bond_level_deformation_gradient( deformation_gradient, ba_deformation_gradient, ) - mean_deformation_gradient = zeros(dof, dof) + mean_deformation_gradient = @MMatrix zeros(dof, dof) for iID in nodes for (jID, nID) in enumerate(nlist[iID]) @views mean_deformation_gradient = diff --git a/test/unit_tests/Core/Solver/ut_Verlet.jl b/test/unit_tests/Core/Solver/ut_Verlet.jl index d87e17c3..6089df85 100644 --- a/test/unit_tests/Core/Solver/ut_Verlet.jl +++ b/test/unit_tests/Core/Solver/ut_Verlet.jl @@ -76,7 +76,7 @@ volume = [0.5, 0.5, 0.5, 0.5, 0.5] density = [1e-6, 1e-6, 3e-6, 3e-6, 1e-6] horizon = [3.1, 3.1, 3.1, 3.1, 3.1] -undeformed_bond = PeriLab.IO.Geometry.bond_geometry( +PeriLab.IO.Geometry.bond_geometry( Vector(1:nnodes), nlist, coor, diff --git a/test/unit_tests/IO/ut_mesh_data.jl b/test/unit_tests/IO/ut_mesh_data.jl index 93ebe4a4..2822a0c6 100644 --- a/test/unit_tests/IO/ut_mesh_data.jl +++ b/test/unit_tests/IO/ut_mesh_data.jl @@ -529,7 +529,7 @@ end test_data_manager.create_constant_bond_field("Bond Geometry", Float64, dof) undeformed_bond_length = test_data_manager.create_constant_bond_field("Bond Length", Float64, 1) - undeformed_bond, undeformed_bond_length = PeriLab.IO.Geometry.bond_geometry( + PeriLab.IO.Geometry.bond_geometry( Vector(1:nnodes), nlist, coor, diff --git a/test/unit_tests/Models/Pre_calculation/ut_pre_bond_associated_correspondence.jl b/test/unit_tests/Models/Pre_calculation/ut_pre_bond_associated_correspondence.jl index a4b23742..4343b229 100644 --- a/test/unit_tests/Models/Pre_calculation/ut_pre_bond_associated_correspondence.jl +++ b/test/unit_tests/Models/Pre_calculation/ut_pre_bond_associated_correspondence.jl @@ -114,7 +114,7 @@ end test_data_manager.create_constant_node_field("Weighted Volume", Float64, 1) - result = Pre_Bond_Associated_Correspondence.compute_weighted_volume( + Pre_Bond_Associated_Correspondence.compute_weighted_volume( nodes, nlist, volume, @@ -124,7 +124,7 @@ end ) expected_weighted_volume = sum(bond_damage[1][:] .* omega[1][:] .* volume[nlist[1]]) - @test isapprox(result[1], expected_weighted_volume) + @test isapprox(weighted_volume[1], expected_weighted_volume) end diff --git a/test/unit_tests/Support/ut_geometry.jl b/test/unit_tests/Support/ut_geometry.jl index 0d248689..42a5a875 100644 --- a/test/unit_tests/Support/ut_geometry.jl +++ b/test/unit_tests/Support/ut_geometry.jl @@ -127,16 +127,22 @@ end coor[4, 1] = 0 coor[4, 2] = 1 - u_bond, u_bond_length = PeriLab.IO.Geometry.calculate_bond_length(1, coor, nlist[1]) + PeriLab.IO.Geometry.calculate_bond_length( + 1, + coor, + nlist[1], + undeformed_bond[1], + undeformed_bond_length[1], + ) - @test u_bond[1, 1] == 0.5 - @test u_bond[1, 2] == 0.5 - @test isapprox(u_bond_length[1], sqrt(0.5)) - @test u_bond[2, 1] == 1 - @test u_bond[2, 2] == 0 - @test u_bond_length[2] == 1 + @test undeformed_bond[1][1, 1] == 0.5 + @test undeformed_bond[1][1, 2] == 0.5 + @test isapprox(undeformed_bond_length[1][1], sqrt(0.5)) + @test undeformed_bond[1][2, 1] == 1 + @test undeformed_bond[1][2, 2] == 0 + @test undeformed_bond_length[1][2] == 1 - undeformed_bond, undeformed_bond_length = PeriLab.IO.Geometry.bond_geometry( + PeriLab.IO.Geometry.bond_geometry( Vector(1:nnodes), nlist, coor, @@ -168,7 +174,7 @@ end @test undeformed_bond[4][1, 1] == 0.5 @test undeformed_bond[4][1, 2] == -0.5 @test undeformed_bond_length[4][1] / sqrt(1.25) - 1 < 1e-8 - undeformed_bond, undeformed_bond_length = PeriLab.IO.Geometry.bond_geometry( + PeriLab.IO.Geometry.bond_geometry( Vector(1:nnodes), nlist, coor, @@ -209,14 +215,13 @@ end undeformed_bond[3][:, :] .= 0 undeformed_bond[4][:, :] .= 0 - undeformed_bond = PeriLab.IO.Geometry.bond_geometry( + PeriLab.IO.Geometry.bond_geometry( Vector(1:nnodes), nlist, coor, undeformed_bond, undeformed_bond_length, ) - @test isnothing(undeformed_bond) end @@ -307,7 +312,7 @@ end coor[4, 1] = 1 coor[4, 2] = 0.5 - undeformed_bond, undeformed_bond_length = PeriLab.IO.Geometry.bond_geometry( + PeriLab.IO.Geometry.bond_geometry( Vector(1:nnodes), nlist, coor, @@ -349,7 +354,7 @@ end deformed_coor[2, 1] = 0.25 deformed_coor[4, 1] = 0.25 - deformed_bond, deformed_bond_length = PeriLab.IO.Geometry.bond_geometry( + PeriLab.IO.Geometry.bond_geometry( Vector(1:nnodes), nlist, deformed_coor, @@ -383,7 +388,7 @@ end deformed_coor[3, 2] = 1.5 deformed_coor[4, 2] = 1.5 - deformed_bond, deformed_bond_length = PeriLab.IO.Geometry.bond_geometry( + PeriLab.IO.Geometry.bond_geometry( Vector(1:nnodes), nlist, deformed_coor, @@ -421,7 +426,7 @@ end deformed_coor[4, 1] = 1.5 deformed_coor[4, 2] = 0.5 - deformed_bond, deformed_bond_length = PeriLab.IO.Geometry.bond_geometry( + PeriLab.IO.Geometry.bond_geometry( Vector(1:nnodes), nlist, deformed_coor, @@ -491,7 +496,7 @@ end inverse_shape_tensor, deformation_gradient, ) - strain = PeriLab.IO.Geometry.compute_strain(nodes, deformation_gradient, strain) + PeriLab.IO.Geometry.compute_strain(nodes, deformation_gradient, strain) for i = 1:nnodes @test strain[i, 1, 1] == 0 @@ -516,7 +521,7 @@ end deformation_gradient_3D[1, 3, 3] = 3.0 strain_3D = test_data_manager.create_constant_node_field("Strain_3D", Float64, "Matrix", 3) - strain_3D = PeriLab.IO.Geometry.compute_strain( + PeriLab.IO.Geometry.compute_strain( view(nodes, eachindex(nodes)), deformation_gradient_3D, strain_3D, @@ -570,7 +575,6 @@ end PeriLab.IO.Geometry.compute_weighted_deformation_gradient( nnodes, - dof, nlist, volume, gradient_weight, From e005734e17319aa5053723fe48f840f952fc048c Mon Sep 17 00:00:00 2001 From: will_cr Date: Wed, 30 Oct 2024 09:52:55 +0100 Subject: [PATCH 3/6] test coverage --- test/unit_tests/Support/ut_helpers.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/unit_tests/Support/ut_helpers.jl b/test/unit_tests/Support/ut_helpers.jl index 4759f9d5..3bc36bbc 100644 --- a/test/unit_tests/Support/ut_helpers.jl +++ b/test/unit_tests/Support/ut_helpers.jl @@ -82,7 +82,10 @@ end end +@testset "ut_get_mapping" begin + @test isnothing(PeriLab.Solver.Helpers.get_mapping(4)) +end @testset "ut_qdim" begin @test isnothing(PeriLab.Solver.Helpers.qdim(0, 2)) From a0d7c9f2dada2eea5e66e00b4d8ed05b4f41cd8c Mon Sep 17 00:00:00 2001 From: will_cr Date: Wed, 30 Oct 2024 10:15:38 +0100 Subject: [PATCH 4/6] memory opt --- src/Models/Pre_calculation/shape_tensor.jl | 4 +- src/Support/geometry.jl | 46 +++++++++++++++++----- src/Support/helpers.jl | 3 +- test/unit_tests/Support/ut_geometry.jl | 6 +-- 4 files changed, 44 insertions(+), 15 deletions(-) diff --git a/src/Models/Pre_calculation/shape_tensor.jl b/src/Models/Pre_calculation/shape_tensor.jl index ed96d7e5..6048952f 100644 --- 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, diff --git a/src/Support/geometry.jl b/src/Support/geometry.jl index 97050d24..b8701032 100644 --- a/src/Support/geometry.jl +++ b/src/Support/geometry.jl @@ -11,7 +11,7 @@ using Rotations include("helpers.jl") using .Helpers: invert, smat, mat_mul_transpose_mat! export bond_geometry -export shape_tensor +export compute_shape_tensors """ @@ -86,7 +86,7 @@ function calculate_bond_length( end end """ - shape_tensor(nodes::Union{SubArray, Vector{Int64}}, nlist, volume, omega, bond_damage, undeformed_bond, shape_tensor, inverse_shape_tensor) + compute_shape_tensors(nodes::Union{SubArray, Vector{Int64}}, nlist, volume, omega, bond_damage, undeformed_bond, shape_tensor, inverse_shape_tensor) Calculate the shape tensor and its inverse for a set of nodes in a computational mechanics context. @@ -122,9 +122,9 @@ undeformed_bond = rand(Float64, length(nodes), length(nlist[1]), dof) shape_tensor = zeros(Float64, length(nodes), dof, dof) inverse_shape_tensor = zeros(Float64, length(nodes), dof, dof) -shape_tensor(nodes, nlist, volume, omega, bond_damage, undeformed_bond, shape_tensor, inverse_shape_tensor) +compute_shape_tensors(nodes, nlist, volume, omega, bond_damage, undeformed_bond, shape_tensor, inverse_shape_tensor) """ -function shape_tensor( +function compute_shape_tensors( nodes::Union{SubArray,Vector{Int64}}, nlist, volume, @@ -137,11 +137,14 @@ function shape_tensor( for iID in nodes - shape_tensor[iID, :, :] = calculate_shape_tensor( - volume[nlist[iID]], + compute_shape_tensor( + shape_tensor, + volume, omega[iID], bond_damage[iID], undeformed_bond[iID], + nlist, + iID, ) #mul!( @@ -155,13 +158,38 @@ function shape_tensor( "Shape Tensor is singular and cannot be inverted).\n - Check if your mesh is 3D, but has only one layer of nodes\n - Check number of damaged bonds.", ) end - return shape_tensor, inverse_shape_tensor + end -function calculate_shape_tensor(volume, omega, bond_damage, undeformed_bond) +#function compute_shape_tensor(volume, omega, bond_damage, undeformed_bond) +# +# return (bond_damage .* volume .* omega .* undeformed_bond)' * undeformed_bond +# +#end - return (bond_damage .* volume .* omega .* undeformed_bond)' * undeformed_bond +function compute_shape_tensor( + shape_tensor, + volume, + omega, + bond_damage, + undeformed_bond, + nlist, + iID, +) + @inbounds @fastmath for m ∈ axes(shape_tensor, 2), n ∈ axes(shape_tensor, 3) + Cmn = zero(eltype(shape_tensor)) + @inbounds @fastmath for k ∈ axes(undeformed_bond, 1) + Cmn += + undeformed_bond[k, m] * + undeformed_bond[k, n] * + volume[nlist[iID][k]] * + omega[k] * + bond_damage[k] + end + shape_tensor[iID, m, n] = Cmn + end + #return (bond_damage .* volume .* omega .* undeformed_bond)' * undeformed_bond end diff --git a/src/Support/helpers.jl b/src/Support/helpers.jl index f9ada13d..ffd332b8 100644 --- a/src/Support/helpers.jl +++ b/src/Support/helpers.jl @@ -45,7 +45,8 @@ function get_mapping(dof::Int64) elseif dof == 3 return (1, 1), (2, 2), (3, 3), (2, 3), (1, 3), (1, 2) else - @error "$dof is no valid mapping option" + @error "$dof is no valid mapping option." + return nothing end end diff --git a/test/unit_tests/Support/ut_geometry.jl b/test/unit_tests/Support/ut_geometry.jl index 42a5a875..4e60ffaf 100644 --- a/test/unit_tests/Support/ut_geometry.jl +++ b/test/unit_tests/Support/ut_geometry.jl @@ -5,8 +5,8 @@ using Test -#include("../../../src/PeriLab.jl") -#using .PeriLab +include("../../../src/PeriLab.jl") +using .PeriLab @testset "ut_compute_bond_level_deformation_gradient" begin nodes = [1] @@ -319,7 +319,7 @@ end undeformed_bond, undeformed_bond_length, ) - shape_tensor, inverse_shape_tensor = PeriLab.IO.Geometry.shape_tensor( + PeriLab.IO.Geometry.compute_shape_tensors( nodes, nlist, volume, From 3ce7b59dd23d3346df36f16579d719be1b40eba6 Mon Sep 17 00:00:00 2001 From: will_cr Date: Thu, 31 Oct 2024 09:16:01 +0100 Subject: [PATCH 5/6] test fix and code improvements --- .../Pre_calculation/deformation_gradient.jl | 7 ++- src/Support/geometry.jl | 58 +++++++++++++----- .../Reference/thermal_flow_correspondence.e | Bin 6876 -> 6876 bytes test/unit_tests/Support/ut_geometry.jl | 15 +++-- 4 files changed, 56 insertions(+), 24 deletions(-) diff --git a/src/Models/Pre_calculation/deformation_gradient.jl b/src/Models/Pre_calculation/deformation_gradient.jl index c2904fe4..ba89d8ea 100644 --- 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,9 +80,10 @@ 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( nodes, + dof, nlist, volume, omega, diff --git a/src/Support/geometry.jl b/src/Support/geometry.jl index b8701032..629a9ed8 100644 --- a/src/Support/geometry.jl +++ b/src/Support/geometry.jl @@ -214,7 +214,7 @@ function bond_associated_deformation_gradient( end """ - compute_deformation_gradient(nodes::Union{SubArray, Vector{Int64}}, nlist, volume, omega, bond_damage, undeformed_bond, deformed_bond, inverse_shape_tensor, deformation_gradient) + compute_deformation_gradients(nodes::Union{SubArray, Vector{Int64}}, nlist, volume, omega, bond_damage, undeformed_bond, deformed_bond, inverse_shape_tensor, deformation_gradient) Calculate the deformation gradient tensor for a set of nodes in a computational mechanics context. @@ -252,8 +252,9 @@ inverse_shape_tensor = rand(Float64, length(nodes), dof, dof) deformation_gradient = zeros(Float64, length(nodes), dof, dof) """ -function compute_deformation_gradient( +function compute_deformation_gradients( nodes::Union{SubArray,Vector{Int64}}, + dof::Int64, nlist::Vector{Vector{Int64}}, volume::Vector{Float64}, omega::Vector{Vector{Float64}}, @@ -263,18 +264,53 @@ function compute_deformation_gradient( inverse_shape_tensor::Array{Float64,3}, deformation_gradient::Array{Float64,3}, ) + if dof == 2 + temp = @MMatrix zeros(2, 2) + else + temp = @MMatrix zeros(3, 3) + end for iID in nodes - deformation_gradient[iID, :, :] = calculate_deformation_gradient( + compute_deformation_gradient( + temp, bond_damage[iID], deformed_bond[iID], undeformed_bond[iID], - volume[nlist[iID]], + volume, omega[iID], + nlist, + iID, ) + dg = @view deformation_gradient[iID, :, :] + mul!(dg, temp, inverse_shape_tensor[iID, :, :]) + end + +end - @views deformation_gradient[iID, :, :] *= inverse_shape_tensor[iID, :, :] +function compute_deformation_gradient( + deformation_gradient, + bond_damage, + deformed_bond, + undeformed_bond, + volume::Union{Vector{Int64},Vector{Float64}}, + omega, + nlist, + iID, +) + #return (bond_damage .* volume .* omega .* deformed_bond)' * undeformed_bond + @inbounds @fastmath for m ∈ axes(deformation_gradient, 1), + n ∈ axes(deformation_gradient, 2) + + Cmn = zero(eltype(deformation_gradient)) + @inbounds @fastmath for k ∈ axes(undeformed_bond, 1) + Cmn += + deformed_bond[k, m] * + undeformed_bond[k, n] * + volume[nlist[iID][k]] * + omega[k] * + bond_damage[k] + end + deformation_gradient[m, n] = Cmn end - return deformation_gradient end @@ -305,16 +341,6 @@ function compute_left_stretch_tensor(deformation_gradient::Matrix{Float64}) return sqrt(deformation_gradient * deformation_gradient') end -function calculate_deformation_gradient( - bond_damage, - deformed_bond, - undeformed_bond, - volume::Union{Vector{Int64},Vector{Float64}}, - omega, -) - return (bond_damage .* volume .* omega .* deformed_bond)' * undeformed_bond - -end function compute_weighted_deformation_gradient( diff --git a/test/fullscale_tests/test_thermal_flow/Reference/thermal_flow_correspondence.e b/test/fullscale_tests/test_thermal_flow/Reference/thermal_flow_correspondence.e index e6b537f16ed17a05ff9178c6fdf1f261d339c6c3..0dfd237bf7a7f4556110b7a1d3376190c00d671e 100644 GIT binary patch delta 1437 zcmca(ddGA_1QXkZ<#k3A>o>$jYP}NqxHOz*3hCtQEo2T=)Fh?G9^6{D%x=<5?KIFQh%ik%u1C4Hn^4kw- zbp%LONBqJ?PqqxH=$#X>bBi}fAH89T1PQaiKISB_` zyLRvzPTJsRVymW|?jamBI6%~5uhdI5K61c8LjpwOdHAG0o_S85+yA`Y^y&dF`jA7j zdv)s7eK={O%&Z5qv|JA1peLUZRhP7I&=vj}rpSgtg9AiIj6YJ%$&=CPV^%Q+eaKFG zIq%9d({a+q*X~`_mfUs-2MrDo!&#?#^POrOaL|wdF=k9@;S-l}@`wugDf>+g7k$XO z|GiH6jcGV(QV9~Q1c+t7tld@hH%{(mEjPR# ze8EK@GG|z~QoH#AZaOELJru=*c!Bit^kq+ZdAO_Bi06fzD9Kuos_x#`6C)ekbp!&yU3uBXH4&oKKo&=^fgP zMsI=gw;a-Q`}ufL)g)Xr&_U{wT@G?w28q+;F=&ueCeIP!k9_Im)xbD2I{||}C{&3CGDz(GR-#F#Oqg-=|@$s;P{r|dT+T=XI9 z{`WfNHzuIblc4-bhfGy=9X%qv1s6RzMqFLe&cP&of1gu61`Q67J>vY4JWlQhucp0P zfk7X#)cx>y)dEvA+8oL^KV%lV^Oe^6N4RKEfXJDB+i5f3o52A)4RnNJhP`FKtld@h zH%{(mEjPR#e8NQ^GG|z~QoH#)8hsJU2husw>{)gqxai4mB-ACZI9U1$hd-Q-L4(7? zN0L9%z{zd9eY@TX4Em5M^URdp3P;iC<52$bLl#+2S1!t*g^LC!5bONS9MfJ$IN+cm z3B-Xf>(4ovI43u!U(a4#{(*}=WNg&-c(MO0oU~ Date: Sun, 3 Nov 2024 01:08:11 +0000 Subject: [PATCH 6/6] CompatHelper: bump compat for Tensorial to 0.18, (keep existing compat) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 483caba0..b68e916d 100644 --- a/Project.toml +++ b/Project.toml @@ -72,7 +72,7 @@ Rotations = "^1" StaticArrays = "^1" Statistics = "^1" TensorOperations = "^4, 5" -Tensorial = "0.16, 0.17" +Tensorial = "0.16, 0.17, 0.18" Tensors = "^1" Test = "^1" TestSetExtensions = "^3"