Added material parameter input through mesh #142 #143

src/Core/Solver/Verlet.jl

@@ -111,7 +111,7 @@ This function depends on the following data fields from the `datamanager` module
 - `get_field("Volume")`: Returns the volume field.
 - `get_field("Horizon")`: Returns the horizon field.
 """
-function compute_mechanical_critical_time_step(nodes::Union{SubArray,Vector{Int64}}, datamanager::Module, bulkModulus::Union{Float64,Int64})
+function compute_mechanical_critical_time_step(nodes::Union{SubArray,Vector{Int64}}, datamanager::Module, bulkModulus::Union{Float64,Int64,SubArray,Vector{Float64}})
     critical_time_step::Float64 = 1.0e50
     nlist = datamanager.get_nlist()
     density = datamanager.get_field("Density")
@@ -122,7 +122,7 @@ function compute_mechanical_critical_time_step(nodes::Union{SubArray,Vector{Int6
     for iID in nodes
         denominator = get_cs_denominator(volume[nlist[iID]], undeformed_bond_length[iID])
 
-        springConstant = 18.0 * bulkModulus / (pi * horizon[iID] * horizon[iID] * horizon[iID] * horizon[iID])
+        springConstant = 18.0 * maximum(bulkModulus) / (pi * horizon[iID] * horizon[iID] * horizon[iID] * horizon[iID])
 
         t = density[iID] / (denominator * springConstant)
         critical_time_step = test_timestep(t, critical_time_step)

src/IO/logging.jl

@@ -243,7 +243,7 @@ function init_logging(filename::String, debug::Bool, silent::Bool, rank::Int64,
         error_logger = FormatLogger(log_file; append=false) do io, args
             if args.level == Logging.Error
                 close_result_files(result_files)
-                throw(ErrorException(args.message))
+                throw(TaskFailedException(args.message))
             end
         end
         filtered_logger = ActiveFilteredLogger(progress_filter, ConsoleLogger(stderr))

src/PeriLab.jl

@@ -276,7 +276,7 @@ function main(filename::String; output_dir::String="", dry_run::Bool=false, verb
         catch e
             if e isa InterruptException
                 @info "PeriLab was interrupted"
-            elseif !isa(e, ErrorException)
+            elseif !isa(e, TaskFailedException)
                 rethrow(e)
             end
         end

src/Physics/Material/Material_Models/Ordinary/Ordinary.jl

@@ -94,7 +94,7 @@ Calculate the symmetry parameters based on the given material symmetry. These pa
 - `gamma::Float64`: Gamma parameter.
 - `kappa::Float64`: Kappa parameter.
 """
-function calculate_symmetry_params(symmetry::String, shear_modulus::Float64, bulk_modulus::Float64)
+function calculate_symmetry_params(symmetry::String, shear_modulus::Union{Float64,SubArray,Vector{Float64}}, bulk_modulus::Union{Float64,SubArray,Vector{Float64}})
     three_bulk_modulus = 3 * bulk_modulus
     # from Peridigm damage model. to be checked with literature
     if symmetry == "plane stress"

src/Physics/Material/Material_Models/PD_Solid_Elastic.jl

@@ -105,7 +105,6 @@ function compute_forces(datamanager::Module, nodes::Union{SubArray,Vector{Int64}
 
     @timeit to "Weighted Volume" weighted_volume = compute_weighted_volume(nodes, nlist, undeformed_bond_length, bond_damage, omega, volume)
     @timeit to "Dilatation" theta = compute_dilatation(nodes, nneighbors, nlist, undeformed_bond_length, deformed_bond_length, bond_damage, volume, weighted_volume, omega)
-
     @timeit to "Bond Forces" bond_force_deviatoric_part, bond_force_isotropic_part = elastic(nodes, dof, undeformed_bond_length, deformed_bond_length, bond_damage, theta, weighted_volume, omega, material_parameter, bond_force_deviatoric_part, bond_force_isotropic_part)
     bond_force = get_bond_forces(nodes, bond_force_deviatoric_part + bond_force_isotropic_part, deformed_bond, deformed_bond_length, bond_force)
 
@@ -139,11 +138,10 @@ function elastic(nodes::Union{SubArray,Vector{Int64}}, dof::Int64, undeformed_bo
 
     shear_modulus = material["Shear Modulus"]
     bulk_modulus = material["Bulk Modulus"]
-
     symmetry::String = get_symmetry(material)
-    kappa::Float64 = 0
-    gamma::Float64 = 0
-    alpha::Float64 = 0
+    # kappa::Float64 = 0
+    # gamma::Float64 = 0
+    # alpha::Float64 = 0
     deviatoric_deformation = @MVector zeros(Float64, dof)
 
     alpha, gamma, kappa = Ordinary.calculate_symmetry_params(symmetry, shear_modulus, bulk_modulus)
@@ -154,8 +152,13 @@ function elastic(nodes::Union{SubArray,Vector{Int64}}, dof::Int64, undeformed_bo
             continue
         end
         deviatoric_deformation = deformed_bond_length[iID] .- undeformed_bond_length[iID] - (gamma * theta[iID] / 3) .* undeformed_bond_length[iID]
-        bond_force_deviatoric_part[iID] = bond_damage[iID] .* omega[iID] .* alpha .* deviatoric_deformation ./ weighted_volume[iID]
-        bond_force_isotropic_part[iID] = bond_damage[iID] .* omega[iID] .* kappa .* theta[iID] .* undeformed_bond_length[iID] ./ weighted_volume[iID]
+        if alpha isa Float64
+            bond_force_deviatoric_part[iID] = bond_damage[iID] .* omega[iID] .* alpha .* deviatoric_deformation ./ weighted_volume[iID]
+            bond_force_isotropic_part[iID] = bond_damage[iID] .* omega[iID] .* kappa .* theta[iID] .* undeformed_bond_length[iID] ./ weighted_volume[iID]
+        else
+            bond_force_deviatoric_part[iID] = bond_damage[iID] .* omega[iID] .* alpha[iID] .* deviatoric_deformation ./ weighted_volume[iID]
+            bond_force_isotropic_part[iID] = bond_damage[iID] .* omega[iID] .* kappa[iID] .* theta[iID] .* undeformed_bond_length[iID] ./ weighted_volume[iID]
+        end
     end
     deviatoric_deformation = nothing
     return bond_force_deviatoric_part, bond_force_isotropic_part

src/Physics/Material/material_basis.jl

@@ -4,14 +4,19 @@
 using LinearAlgebra
 using StaticArrays
 
-function get_value(datamanager::Module, parameter::Union{Dict{Any,Any},Dict{String,Any}}, key::String, field_allocated::Bool)
+function get_value(datamanager::Module, parameter::Union{Dict{Any,Any},Dict{String,Any}}, any_field_allocated::Bool, key::String, field_allocated::Bool)
     if field_allocated
-        return datamanager.get_field(key)
+        return datamanager.get_field(replace(key, " " => "_"))
     end
     if haskey(parameter, key)
-        return parameter[key]
+        if any_field_allocated
+            return datamanager.create_constant_node_field(replace(key, " " => "_"), Float64, 1, parameter[key])
+        else
+            return parameter[key]
+        end
     end
-    return nothing
+
+    return datamanager.create_constant_node_field(replace(key, " " => "_"), Float64, 1)
 end
 
 """
@@ -28,26 +33,23 @@ function get_all_elastic_moduli(datamanager::Module, parameter::Union{Dict{Any,A
             return nothing
         end
     end
-    bulk = false
-    Youngs = false
-    shear = false
-    Poissons = false
-    field_allocated = false
 
-    bulk_field = datamanager.has_key("Bulk Modulus")
-    Youngs_field = datamanager.has_key("Young's Modulus")
-    Poissons_field = datamanager.has_key("Poisson's Ratio")
-    shear_field = datamanager.has_key("Shear Modulus")
+    bulk_field = datamanager.has_key("Bulk_Modulus")
+    Youngs_field = datamanager.has_key("Young's_Modulus")
+    Poissons_field = datamanager.has_key("Poisson's_Ratio")
+    shear_field = datamanager.has_key("Shear_Modulus")
+
+    any_field_allocated = bulk_field | Youngs_field | Poissons_field | shear_field
 
     bulk = haskey(parameter, "Bulk Modulus") | bulk_field
     Youngs = haskey(parameter, "Young's Modulus") | Youngs_field
     Poissons = haskey(parameter, "Poisson's Ratio") | Poissons_field
     shear = haskey(parameter, "Shear Modulus") | shear_field
 
-    K = get_value(datamanager, parameter, "Bulk Modulus", bulk_field)
-    E = get_value(datamanager, parameter, "Young's Modulus", Youngs_field)
-    G = get_value(datamanager, parameter, "Shear Modulus", shear_field)
-    nu = get_value(datamanager, parameter, "Poisson's Ratio", Poissons_field)
+    K = get_value(datamanager, parameter, any_field_allocated, "Bulk Modulus", bulk_field)
+    E = get_value(datamanager, parameter, any_field_allocated, "Young's Modulus", Youngs_field)
+    G = get_value(datamanager, parameter, any_field_allocated, "Shear Modulus", shear_field)
+    nu = get_value(datamanager, parameter, any_field_allocated, "Poisson's Ratio", Poissons_field)
 
     if bulk && Poissons
         E = 3 .* K .* (1 .- 2 .* nu)
@@ -78,12 +80,17 @@ function get_all_elastic_moduli(datamanager::Module, parameter::Union{Dict{Any,A
     if bulk + Youngs + shear + Poissons < 2
         @error "Minimum of two parameters are needed for isotropic material"
     end
-
     parameter["Bulk Modulus"] = K
     parameter["Young's Modulus"] = E
     parameter["Shear Modulus"] = G
     parameter["Poisson's Ratio"] = nu
     parameter["Computed"] = true
+    if any_field_allocated
+        datamanager.get_field("Bulk_Modulus") .= K
+        datamanager.get_field("Young's_Modulus") .= E
+        datamanager.get_field("Shear_Modulus") .= G
+        datamanager.get_field("Poisson's_Ratio") .= nu
+    end
 end
 
 """

test/fullscale_tests/test_material_field/Reference/material_field.e.license

@@ -0,0 +1,3 @@
+SPDX-FileCopyrightText: 2023 Christian Willberg <christian.willberg@dlr.de>, Jan-Timo Hesse <jan-timo.hesse@dlr.de>
+
+SPDX-License-Identifier: BSD-3-Clause

test/fullscale_tests/test_material_field/material_field.yaml

@@ -0,0 +1,57 @@
+# SPDX-FileCopyrightText: 2023 Christian Willberg <christian.willberg@dlr.de>, Jan-Timo Hesse <jan-timo.hesse@dlr.de>
+#
+# 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
+    Type: "Text File"
+    Input Mesh File: "test_mesh.txt"
+  Physics:
+    Material Models:
+      Mat_1:
+        Material Model: "PD Solid Elastic"
+        Bulk Modulus: 2.5e+3
+        Shear Modulus: 1.15e3
+      Mat_2:
+        Material Model: "PD Solid Elastic"
+        Bulk Modulus: 3.5e+3
+        Shear Modulus: 2.15e3
+  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*y"
+    BC_2:
+      Variable: "Displacements"
+      Node Set: "Node Set 1"
+      Coordinate: "y"
+      Value: "0.1*x"
+  Solver:
+    Initial Time: 0.0
+    Final Time: 5.0e-6
+    Verlet:
+      Safety Factor: 1.0
+  Outputs:
+    Output1:
+      Output Filename: "material_field"
+      Output Frequency: 1
+      Output Variables:
+        Displacements: True
+        Forces: True
+        Volume: True
+        Bulk_Modulus: True
+        Shear_Modulus: True
+        Young's_Modulus: True
+        Poisson's_Ratio: True

test/fullscale_tests/test_material_field/test_material_field.cmd

@@ -0,0 +1,7 @@
+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

test/fullscale_tests/test_material_field/test_material_field.cmd.license

@@ -0,0 +1,3 @@
+SPDX-FileCopyrightText: 2023 Christian Willberg <christian.willberg@dlr.de>, Jan-Timo Hesse <jan-timo.hesse@dlr.de>
+
+SPDX-License-Identifier: BSD-3-Clause

test/fullscale_tests/test_material_field/test_material_field.jl

@@ -0,0 +1,9 @@
+# SPDX-FileCopyrightText: 2023 Christian Willberg <christian.willberg@dlr.de>, Jan-Timo Hesse <jan-timo.hesse@dlr.de>
+#
+# SPDX-License-Identifier: BSD-3-Clause
+
+
+folder_name = basename(@__FILE__)[1:end-3]
+cd("fullscale_tests/" * folder_name) do
+    run_perilab("material_field", 1, true, folder_name)
+end

test/fullscale_tests/test_material_field/test_mesh.txt

@@ -0,0 +1,17 @@
+# SPDX-FileCopyrightText: 2023 Christian Willberg <christian.willberg@dlr.de>, Jan-Timo Hesse <jan-timo.hesse@dlr.de>
+#
+# SPDX-License-Identifier: BSD-3-Clause
+
+header: x y block_id volume Bulk_Modulus
+-1.5 0 1 1.e-01 2.5e+3
+-0.5 0 1 1.e-01 2.5e+3
+0.5 0 1 1.e-01 3.5e+3
+1.5 0 1 1.e-01 3.5e+3
+-1.5 1 1 1.e-01 4.5e+3
+-0.5 1 1 1.e-01 4.5e+3
+0.5 1 2 1.e-01 2.5e+3
+1.5 1 2 1.e-01 2.5e+3
+-1.5 2 1 1.e-01 3.5e+3
+-0.5 2 1 1.e-01 3.5e+3
+0.5 2 2 1.e-01 4.5e+3
+1.5 2 2 1.e-01 4.5e+3

test/runtests.jl

@@ -229,6 +229,9 @@ MPI.Init()
         @testset "test_aniso_damage" begin
             @includetests["fullscale_tests/test_aniso_damage/test_aniso_damage"]
         end
+        @testset "test_material_field" begin
+            @includetests["fullscale_tests/test_material_field/test_material_field"]
+        end
         @testset "test_FEM" begin
             @includetests["fullscale_tests/test_FEM/test_FEM"]
         end