diff --git a/src/Core/Solver/Verlet.jl b/src/Core/Solver/Verlet.jl index 790674c6..90ee626f 100644 --- a/src/Core/Solver/Verlet.jl +++ b/src/Core/Solver/Verlet.jl @@ -220,6 +220,12 @@ function compute_crititical_time_step( if mechanical bulk_modulus = datamanager.get_property(iblock, "Material Model", "Bulk Modulus") + g_xy = datamanager.get_property(iblock, "Material Model", "Shear Modulus XY") + g_yz = datamanager.get_property(iblock, "Material Model", "Shear Modulus YZ") + g_zx = datamanager.get_property(iblock, "Material Model", "Shear Modulus ZX") + c_44 = datamanager.get_property(iblock, "Material Model", "C44") + c_55 = datamanager.get_property(iblock, "Material Model", "C55") + c_66 = datamanager.get_property(iblock, "Material Model", "C66") if !isnothing(bulk_modulus) t = compute_mechanical_critical_time_step( block_nodes[iblock], @@ -227,6 +233,21 @@ function compute_crititical_time_step( bulk_modulus, ) critical_time_step = test_timestep(t, critical_time_step) + elseif !isnothing(g_xy) && !isnothing(g_yz) && !isnothing(g_zx) + t = compute_mechanical_critical_time_step( + block_nodes[iblock], + datamanager, + maximum([g_xy, g_yz, g_zx]), + ) + critical_time_step = test_timestep(t, critical_time_step) + #TODO: temporary solution!!! + elseif !isnothing(c_44) && !isnothing(c_55) && !isnothing(c_66) + t = compute_mechanical_critical_time_step( + block_nodes[iblock], + datamanager, + maximum([c_44 / 2, c_55 / 2, c_66 / 2]), + ) + critical_time_step = test_timestep(t, critical_time_step) else @error "No time step for material is determined because of missing properties." return nothing diff --git a/src/Models/Material/Material_Basis.jl b/src/Models/Material/Material_Basis.jl index 1270dbe5..98a4fdfb 100644 --- a/src/Models/Material/Material_Basis.jl +++ b/src/Models/Material/Material_Basis.jl @@ -131,6 +131,36 @@ function get_all_elastic_moduli( nu = nu_fixed poissons = true end + if haskey(parameter, "Symmetry") && + occursin("anisotropic", lowercase(parameter["Symmetry"])) + E_x = haskey(parameter, "Young's Modulus X") + E_y = haskey(parameter, "Young's Modulus Y") + E_z = haskey(parameter, "Young's Modulus Z") + nu_xy = haskey(parameter, "Poisson's Ratio XY") + nu_yz = haskey(parameter, "Poisson's Ratio YZ") + nu_zx = haskey(parameter, "Poisson's Ratio ZX") + g_xy = haskey(parameter, "Shear Modulus XY") + g_yz = haskey(parameter, "Shear Modulus YZ") + g_zx = haskey(parameter, "Shear Modulus ZX") + if E_x + E_y + E_z + nu_xy + nu_yz + nu_zx + g_xy + g_yz + g_zx > 0 + if !E_x || !E_y || !E_z || !nu_xy || !nu_yz || !nu_zx || !g_xy || !g_yz || !g_zx + @error "Non isotropic material requires Young's Modulus X, Y, Z, Poisson's Ratio XY, YZ, ZX, Shear Modulus XY, YZ, ZX" + return nothing + end + parameter["Compute_Hook"] = true + else + for iID = 1:6 + for jID = iID:6 + if !("C" * string(iID) * string(jID) in keys(parameter)) + @error "C" * string(iID) * string(jID) * " not defined" + return nothing + end + end + end + parameter["Compute_Hook"] = false + end + return + end # tbd non isotropic material check if bulk + youngs + shear + poissons < 2 @@ -193,21 +223,59 @@ Returns the Hooke matrix of the material. """ function get_Hooke_matrix(parameter::Dict, symmetry::String, dof::Int64, ID::Int64 = 1) """https://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_plane_stress.cfm""" - - if occursin("anisotropic", symmetry) + if occursin("anisotropic", lowercase(symmetry)) aniso_matrix = get_MMatrix(36) - for iID = 1:6 - for jID = iID:6 - if "C" * string(iID) * string(jID) in keys(parameter) + + if parameter["Compute_Hook"] + + E_x = parameter["Young's Modulus X"] + E_y = parameter["Young's Modulus Y"] + E_z = parameter["Young's Modulus Z"] + nu_xy = parameter["Poisson's Ratio XY"] + nu_yz = parameter["Poisson's Ratio YZ"] + nu_zx = parameter["Poisson's Ratio ZX"] + g_xy = parameter["Shear Modulus XY"] + g_yz = parameter["Shear Modulus YZ"] + g_zx = parameter["Shear Modulus ZX"] + + nu_yx = nu_xy * E_y / E_x + nu_xz = nu_zx * E_x / E_z + nu_zy = nu_yz * E_z / E_y + + delta = + ( + 1 - nu_xy * nu_yx - nu_yz * nu_zy - nu_zx * nu_xz - + 2 * nu_xy * nu_yz * nu_zx + ) / E_x * + E_y * + E_z + + aniso_matrix[1, 1] = (1 - nu_yz * nu_zy) / E_y * E_z * delta + aniso_matrix[2, 2] = (1 - nu_zx * nu_xz) / E_z * E_x * delta + aniso_matrix[3, 3] = (1 - nu_xy * nu_yx) / E_x * E_y * delta + + aniso_matrix[1, 2] = (nu_yx + nu_zx * nu_yz) / E_y * E_z * delta + aniso_matrix[2, 1] = (nu_xy + nu_xz * nu_zy) / E_z * E_x * delta + + aniso_matrix[1, 3] = (nu_zx + nu_yx * nu_zy) / E_y * E_z * delta + aniso_matrix[3, 1] = (nu_xz + nu_xy * nu_yz) / E_x * E_y * delta + + aniso_matrix[2, 3] = (nu_zy + nu_zx * nu_xy) / E_z * E_x * delta + aniso_matrix[3, 2] = (nu_yz + nu_xz * nu_yx) / E_x * E_y * delta + + aniso_matrix[4, 4] = 2 * g_yz + aniso_matrix[5, 5] = 2 * g_zx + aniso_matrix[6, 6] = 2 * g_xy + else + for iID = 1:6 + for jID = iID:6 value = parameter["C"*string(iID)*string(jID)] - else - @error "C" * string(iID) * string(jID) * " not defined" - return nothing + aniso_matrix[iID, jID] = value + aniso_matrix[jID, iID] = value end - aniso_matrix[iID, jID] = value - aniso_matrix[jID, iID] = value end end + if dof == 3 return aniso_matrix elseif occursin("plane strain", symmetry) @@ -218,12 +286,16 @@ function get_Hooke_matrix(parameter::Dict, symmetry::String, dof::Int64, ID::Int matrix[3, 3] = aniso_matrix[6, 6] return matrix elseif occursin("plane stress", symmetry) - matrix = get_MMatrix(36) + @info aniso_matrix inv_aniso = invert(aniso_matrix, "Hooke matrix not invertable") + @info inv_aniso + matrix = get_MMatrix(36) + @info aniso_matrix matrix[1:2, 1:2] = inv_aniso[1:2, 1:2] matrix[3, 1:2] = inv_aniso[6, 1:2] matrix[1:2, 3] = inv_aniso[1:2, 6] matrix[3, 3] = inv_aniso[6, 6] + @info matrix return invert(matrix, "Hooke matrix not invertable") else @error "2D model defintion is missing; plane stress or plane strain " diff --git a/src/Support/Helpers.jl b/src/Support/Helpers.jl index b73a66d5..c7ad908f 100644 --- a/src/Support/Helpers.jl +++ b/src/Support/Helpers.jl @@ -112,15 +112,18 @@ function add_in_place!( return C end function get_MMatrix(len::Int64) + global A2x2 + global A3x3 + global A6x6 if len == 4 - global A2x2 + A2x2 .= 0 return A2x2 elseif len == 9 - global A3x3 + A3x3 .= 0 return A3x3 elseif len == 36 - global A6x6 + A6x6 .= 0 return A6x6 else @error "MMatrix length $len not pre-allocated. Please add your size to helper.jl in get_MMatrix." diff --git a/src/Support/Parameters/parameter_handling.jl b/src/Support/Parameters/parameter_handling.jl index 0038ad00..6010390b 100644 --- a/src/Support/Parameters/parameter_handling.jl +++ b/src/Support/Parameters/parameter_handling.jl @@ -202,9 +202,24 @@ global expected_structure = Dict( "Bond Associated" => [Bool, false], "Bond Horizon" => [Union{Float64,Int64}, false], "Poisson's Ratio" => [Union{Float64,Int64}, false], + "Poisson's Ratio XY" => + [Union{Float64,Int64}, false], + "Poisson's Ratio YZ" => + [Union{Float64,Int64}, false], + "Poisson's Ratio ZX" => + [Union{Float64,Int64}, false], "Young's Modulus" => [Union{Float64,Int64}, false], + "Young's Modulus X" => + [Union{Float64,Int64}, false], + "Young's Modulus Y" => + [Union{Float64,Int64}, false], + "Young's Modulus Z" => + [Union{Float64,Int64}, false], "Bulk Modulus" => [Union{Float64,Int64}, false], "Shear Modulus" => [Union{Float64,Int64}, false], + "Shear Modulus XY" => [Union{Float64,Int64}, false], + "Shear Modulus YZ" => [Union{Float64,Int64}, false], + "Shear Modulus ZX" => [Union{Float64,Int64}, false], "Yield Stress" => [Union{Float64,Int64}, false], "Zero Energy Control" => [String, false], "C11" => [Union{Float64,Int64}, false], diff --git a/test/fullscale_tests/test_aniso/Reference/aniso_eng.e b/test/fullscale_tests/test_aniso/Reference/aniso_eng.e new file mode 100644 index 00000000..1f45d3f5 Binary files /dev/null and b/test/fullscale_tests/test_aniso/Reference/aniso_eng.e differ diff --git a/test/fullscale_tests/test_aniso/Reference/aniso_eng.e.license b/test/fullscale_tests/test_aniso/Reference/aniso_eng.e.license new file mode 100644 index 00000000..2b47ef51 --- /dev/null +++ b/test/fullscale_tests/test_aniso/Reference/aniso_eng.e.license @@ -0,0 +1,3 @@ +SPDX-FileCopyrightText: 2023 Christian Willberg , Jan-Timo Hesse + +SPDX-License-Identifier: BSD-3-Clause diff --git a/test/fullscale_tests/test_aniso/Reference/aniso_hooke.e b/test/fullscale_tests/test_aniso/Reference/aniso_hooke.e new file mode 100644 index 00000000..0217e424 Binary files /dev/null and b/test/fullscale_tests/test_aniso/Reference/aniso_hooke.e differ diff --git a/test/fullscale_tests/test_aniso/Reference/aniso_hooke.e.license b/test/fullscale_tests/test_aniso/Reference/aniso_hooke.e.license new file mode 100644 index 00000000..2b47ef51 --- /dev/null +++ b/test/fullscale_tests/test_aniso/Reference/aniso_hooke.e.license @@ -0,0 +1,3 @@ +SPDX-FileCopyrightText: 2023 Christian Willberg , Jan-Timo Hesse + +SPDX-License-Identifier: BSD-3-Clause diff --git a/test/fullscale_tests/test_aniso/aniso_eng.yaml b/test/fullscale_tests/test_aniso/aniso_eng.yaml new file mode 100644 index 00000000..37fdeead --- /dev/null +++ b/test/fullscale_tests/test_aniso/aniso_eng.yaml @@ -0,0 +1,80 @@ +# SPDX-FileCopyrightText: 2023 Christian Willberg , Jan-Timo Hesse +# +# SPDX-License-Identifier: BSD-3-Clause + +PeriLab: + Discretization: + Node Sets: + Node Set 1: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 + Type: "Text File" + Input Mesh File: "test_mesh.txt" + Models: + Material Models: + Mat_1: + Material Model: "Correspondence Elastic" + Symmetry: Anisotropic + Young's Modulus X: 2.5e3 + Young's Modulus Y: 2.5e3 + Young's Modulus Z: 2.5e3 + Poisson's Ratio XY: 0.3 + Poisson's Ratio YZ: 0.3 + Poisson's Ratio ZX: 0.3 + Shear Modulus XY: 2.15e3 + Shear Modulus YZ: 2.15e3 + Shear Modulus ZX: 2.15e3 + Zero Energy Control: "Global" + Mat_2: + Material Model: "Correspondence Elastic" + Symmetry: Anisotropic + Young's Modulus X: 2.5e3 + Young's Modulus Y: 2.5e3 + Young's Modulus Z: 2.5e3 + Poisson's Ratio XY: 0.3 + Poisson's Ratio YZ: 0.3 + Poisson's Ratio ZX: 0.3 + Shear Modulus XY: 2.15e3 + Shear Modulus YZ: 2.15e3 + Shear Modulus ZX: 2.15e3 + Zero Energy Control: "Global" + Blocks: + block_1: + Material Model: "Mat_1" + Density: 2.7e-9 + Horizon: 2.1 + block_2: + Material Model: "Mat_2" + Density: 2.7e-9 + Horizon: 2.1 + Boundary Conditions: + BC_1: + Variable: "Displacements" + Node Set: "Node Set 1" + Coordinate: "x" + Value: "0.1*x" + Type: Dirichlet + BC_2: + Variable: "Displacements" + Node Set: "Node Set 1" + Coordinate: "y" + Value: "0.0" + Type: Dirichlet + BC_3: + Variable: "Displacements" + Node Set: "Node Set 1" + Coordinate: "z" + Value: "0.0" + Type: Dirichlet + Solver: + Initial Time: 0.0 + Final Time: 5.0e-6 + Verlet: + Safety Factor: 1.0 + Outputs: + Output1: + Output Filename: "aniso_eng" + Output File Type: Exodus + Output Frequency: 1 + Output Variables: + Displacements: True + Forces: True + Volume: True diff --git a/test/fullscale_tests/test_aniso/aniso_hooke.yaml b/test/fullscale_tests/test_aniso/aniso_hooke.yaml new file mode 100644 index 00000000..b2ec682e --- /dev/null +++ b/test/fullscale_tests/test_aniso/aniso_hooke.yaml @@ -0,0 +1,104 @@ +# SPDX-FileCopyrightText: 2023 Christian Willberg , Jan-Timo Hesse +# +# SPDX-License-Identifier: BSD-3-Clause + +PeriLab: + Discretization: + Node Sets: + Node Set 1: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 + Type: "Text File" + Input Mesh File: "test_mesh.txt" + Models: + Material Models: + Mat_1: + Material Model: "Correspondence Elastic" + Symmetry: Anisotropic + C11: 53619.92831989742 + C12: 42419.92831989742 + C13: 5826.105142734548 + C14: 0.0 + C15: 0.0 + C16: 39393.823177162834 + C22: 53619.92831989742 + C23: 5826.105142734548 + C24: 0.0 + C25: 0.0 + C26: 39393.823177162834 + C33: 11402.506122303113 + C34: 0.0 + C35: 0.0 + C36: 23.599020431513775 + C44: 4200.000000000001 + C45: 1400.0 + C46: 0.0 + C55: 4200.000000000001 + C56: 0.0 + C66: 42170.22415673136 + Zero Energy Control: "Global" + Mat_2: + Material Model: "Correspondence Elastic" + Symmetry: Anisotropic + C11: 53619.92831989742 + C12: 42419.92831989742 + C13: 5826.105142734548 + C14: 0.0 + C15: 0.0 + C16: 39393.823177162834 + C22: 53619.92831989742 + C23: 5826.105142734548 + C24: 0.0 + C25: 0.0 + C26: 39393.823177162834 + C33: 11402.506122303113 + C34: 0.0 + C35: 0.0 + C36: 23.599020431513775 + C44: 4200.000000000001 + C45: 1400.0 + C46: 0.0 + C55: 4200.000000000001 + C56: 0.0 + C66: 42170.22415673136 + Zero Energy Control: "Global" + Blocks: + block_1: + Material Model: "Mat_1" + Density: 2.7e-9 + Horizon: 2.1 + block_2: + Material Model: "Mat_2" + Density: 2.7e-9 + Horizon: 2.1 + Boundary Conditions: + BC_1: + Variable: "Displacements" + Node Set: "Node Set 1" + Coordinate: "x" + Value: "0.1*x" + Type: Dirichlet + BC_2: + Variable: "Displacements" + Node Set: "Node Set 1" + Coordinate: "y" + Value: "0.0" + Type: Dirichlet + BC_3: + Variable: "Displacements" + Node Set: "Node Set 1" + Coordinate: "z" + Value: "0.0" + Type: Dirichlet + Solver: + Initial Time: 0.0 + Final Time: 5.0e-6 + Verlet: + Safety Factor: 1.0 + Outputs: + Output1: + Output Filename: "aniso_hooke" + Output File Type: Exodus + Output Frequency: 1 + Output Variables: + Displacements: True + Forces: True + Volume: True diff --git a/test/fullscale_tests/test_aniso/test_aniso.cmd b/test/fullscale_tests/test_aniso/test_aniso.cmd new file mode 100755 index 00000000..56184424 --- /dev/null +++ b/test/fullscale_tests/test_aniso/test_aniso.cmd @@ -0,0 +1,10 @@ +DEFAULT TOLERANCE absolute 1.0E-9 +COORDINATES absolute 1.0E-12 +TIME STEPS absolute 1.0E-14 +NODAL VARIABLES absolute 1.0E-12 + DisplacementsX absolute 1.0E-9 + DisplacementsY absolute 1.0E-9 + DisplacementsZ absolute 1.0E-9 + ForcesX absolute 5.0E-8 + ForcesY absolute 5.0E-8 + ForcesZ absolute 5.0E-8 diff --git a/test/fullscale_tests/test_aniso/test_aniso.cmd.license b/test/fullscale_tests/test_aniso/test_aniso.cmd.license new file mode 100644 index 00000000..2b47ef51 --- /dev/null +++ b/test/fullscale_tests/test_aniso/test_aniso.cmd.license @@ -0,0 +1,3 @@ +SPDX-FileCopyrightText: 2023 Christian Willberg , Jan-Timo Hesse + +SPDX-License-Identifier: BSD-3-Clause diff --git a/test/fullscale_tests/test_aniso/test_aniso.jl b/test/fullscale_tests/test_aniso/test_aniso.jl new file mode 100644 index 00000000..f657f5da --- /dev/null +++ b/test/fullscale_tests/test_aniso/test_aniso.jl @@ -0,0 +1,11 @@ +# SPDX-FileCopyrightText: 2023 Christian Willberg , Jan-Timo Hesse +# +# SPDX-License-Identifier: BSD-3-Clause + + + +folder_name = basename(@__FILE__)[1:end-3] +cd("fullscale_tests/" * folder_name) do + run_perilab("aniso_eng", 1, true, folder_name) + run_perilab("aniso_hooke", 1, true, folder_name) +end diff --git a/test/fullscale_tests/test_aniso/test_mesh.txt b/test/fullscale_tests/test_aniso/test_mesh.txt new file mode 100644 index 00000000..d24df073 --- /dev/null +++ b/test/fullscale_tests/test_aniso/test_mesh.txt @@ -0,0 +1,41 @@ +# SPDX-FileCopyrightText: 2023 Christian Willberg , Jan-Timo Hesse +# +# SPDX-License-Identifier: BSD-3-Clause + +header: x y z block_id volume +-1.5 0 0 1 1.e-01 +-0.5 0 0 1 1.e-01 +0.5 0 0 1 1.e-01 +1.5 0 0 1 1.e-01 +-1.5 1 0 1 1.e-01 +-0.5 1 0 1 1.e-01 +0.5 1 0 2 1.e-01 +1.5 1 0 2 1.e-01 +-1.5 2 0 1 1.e-01 +-0.5 2 0 1 1.e-01 +0.5 2 0 2 1.e-01 +1.5 2 0 2 1.e-01 +-1.5 0 0.5 1 1.e-01 +-0.5 0 0.5 1 1.e-01 +0.5 0 0.5 1 1.e-01 +1.5 0 0.5 1 1.e-01 +-1.5 1 0.5 1 1.e-01 +-0.5 1 0.5 1 1.e-01 +0.5 1 0.5 2 1.e-01 +1.5 1 0.5 2 1.e-01 +-1.5 2 0.5 1 1.e-01 +-0.5 2 0.5 1 1.e-01 +0.5 2 0.5 2 1.e-01 +1.5 2 0.5 2 1.e-01 +-1.5 0 1 1 1.e-01 +-0.5 0 1 1 1.e-01 +0.5 0 1 1 1.e-01 +1.5 0 1 1 1.e-01 +-1.5 1 1 1 1.e-01 +-0.5 1 1 1 1.e-01 +0.5 1 1 2 1.e-01 +1.5 1 1 2 1.e-01 +-1.5 2 1 1 1.e-01 +-0.5 2 1 1 1.e-01 +0.5 2 1 2 1.e-01 +1.5 2 1 2 1.e-01 diff --git a/test/runtests.jl b/test/runtests.jl index d6e5d99c..8d07b072 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -215,6 +215,9 @@ MPI.Init() @testset "test_two_block" begin include("fullscale_tests/test_two_block/test_two_block.jl") end + @testset "test_aniso" begin + include("fullscale_tests/test_aniso/test_aniso.jl") + end @testset "test_additive_simple" begin include("fullscale_tests/test_additive/test_additive.jl") end diff --git a/test/unit_tests/Models/Material/ut_material_basis.jl b/test/unit_tests/Models/Material/ut_material_basis.jl index 7235f47e..dc18b042 100644 --- a/test/unit_tests/Models/Material/ut_material_basis.jl +++ b/test/unit_tests/Models/Material/ut_material_basis.jl @@ -341,6 +341,7 @@ end "Bulk Modulus" => 5, "Shear Modulus" => 1.25, "Poisson's Ratio" => 0.2, + "Compute_Hook" => false, ) get_all_elastic_moduli(test_data_manager, parameter) @@ -417,6 +418,7 @@ end end end end + symmetry = "anisotropic plane strain" C = get_Hooke_matrix(parameter, symmetry, 2) for iID = 1:2 @@ -433,28 +435,46 @@ end @test C[3, 1] == parameter["C16"] @test C[3, 2] == parameter["C26"] - symmetry = "anisotropic plane stress" - C = get_Hooke_matrix(parameter, symmetry, 2) - for iID = 1:3 - for jID = 1:3 - if C2D_test[iID, jID] != 0 - @test C[iID, jID] / C2D_test[iID, jID] - 1 < 1e-7 - end - end - end + # symmetry = "anisotropic plane stress" + # C = get_Hooke_matrix(parameter, symmetry, 2) + # for iID = 1:3 + # for jID = 1:3 + # if C2D_test[iID, jID] != 0 + # @test C[iID, jID] / C2D_test[iID, jID] - 1 < 1e-7 + # end + # end + # end + #TODO: Check above + + symmetry = "anisotropic" parameter = Dict{String,Any}( - "C11" => 2.0, - "C12" => 3.0, - "C13" => 4.0, - "C22" => 5.0, - "C33" => 7.0, + "Material Model" => "PD Solid Elastic", + "Young's Modulus X" => 5, + "Young's Modulus Y" => 6, + "Young's Modulus Z" => 7, + "Poisson's Ratio XY" => 0.1, + "Poisson's Ratio YZ" => 0.2, + "Poisson's Ratio ZX" => 0.3, + "Shear Modulus XY" => 1, + "Shear Modulus YZ" => 2, + "Shear Modulus ZX" => 3, + "Compute_Hook" => true, ) - @test isnothing(get_Hooke_matrix(parameter, symmetry, 2)) - - symmetry = "anisotropic missing" - @test isnothing(get_Hooke_matrix(parameter, symmetry, 2)) - + C = get_Hooke_matrix(parameter, symmetry, 3) + @info C + @test C[1, 1] == 8.081851555555554 + @test C[1, 2] == 1.525944 + @test C[1, 3] == 2.7806090666666665 + @test C[2, 1] == 0.7785428571428572 + @test C[2, 2] == 4.856624489795917 + @test C[2, 3] == 1.3667752380952378 + @test C[3, 1] == 2.0428964571428567 + @test C[3, 2] == 1.9681563428571425 + @test C[3, 3] == 8.615043839999998 + @test C[4, 4] == 4 + @test C[5, 5] == 6 + @test C[6, 6] == 2 end @testset "ut_compute_Piola_Kirchhoff_stress" begin