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

src/IO/mesh_data.jl

@@ -283,7 +283,7 @@ function get_bond_geometry(datamanager::Module)
     coor = datamanager.get_field("Coordinates")
     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), dof, nlist, coor, undeformed_bond, undeformed_bond_length)
+    undeformed_bond, undeformed_bond_length = Geometry.bond_geometry(Vector(1:nnodes), nlist, coor, undeformed_bond, undeformed_bond_length)
     return datamanager
 end
 

src/Physics/Material/Material_Factory.jl

@@ -72,8 +72,7 @@ function init_material_model(datamanager::Module, nodes::Union{SubArray,Vector{I
 
     if occursin("Correspondence", model_param["Material Model"])
         datamanager.set_model_module("Correspondence", Correspondence)
-        datamanager = Correspondence.init_material_model(datamanager, nodes, model_param)
-        return datamanager
+        return Correspondence.init_material_model(datamanager, nodes, model_param)
     end
 
     material_models = split(model_param["Material Model"], "+")

src/Physics/Material/Material_Models/Bond_Associated_Correspondence.jl

@@ -20,16 +20,16 @@ function init_material_model(datamanager::Module, nodes::Union{SubArray,Vector{I
   bond_horizon = material_parameter["Bond Horizon"]
   weighted_volume = datamanager.create_constant_bond_field("Bond Weighted Volume", Float64, 1)
 
-  try
-    for iID in nodes
-      for (jID, nID) in enumerate(nlist)
-        find_local_neighbors(coordinates[nID, :], coordinates[nlist, :], nlist, bond_horizon)
+  for iID in nodes
+    for nID in nlist[iID]
+      if length(find_local_neighbors(nID, coordinates, nlist[iID], bond_horizon)) == 0
+        @error "Bond horizon is to small. No neighbors were found."
+        return nothing
       end
+
     end
-  catch
-    @error "Bond horizon is to small. Now neighbors were found."
-    return nothing
   end
+  return datamanager
 end
 
 function compute_stress_integral(nodes::Union{SubArray,Vector{Int64}}, nlist::SubArray, omega::SubArray, bond_damage::SubArray, volume::SubArray, weighted_volume::SubArray, bond_geometry::SubArray, bond_length::SubArray, bond_stresses::SubArray, stress_integral::SubArray)
@@ -85,12 +85,12 @@ https://link.springer.com/article/10.1007/s10409-021-01055-5
 
 
 
-function find_local_neighbors(neighbor_coordinate::Union{Vector{Float64},Vector{Int64}}, coordinates::Union{Matrix{Float64},Matrix{Int64}}, nlist::Union{Vector{Int64},SubArray{Int64}}, bond_horizon::Union{Float64,Int64})
+function find_local_neighbors(nID::Int64, coordinates::Union{SubArray,Matrix{Float64},Matrix{Int64}}, nlist::Union{Vector{Int64},SubArray{Int64}}, bond_horizon::Union{Float64,Int64})
   # excludes right now iID node in the coordinates list. Because it is an abritrary sublist it should be fine.
   # saving faster than recalculation?
-  balltree = BallTree(transpose(coordinates))
-
-  return nlist[inrange(balltree, neighbor_coordinate, bond_horizon, true)]
+  nlist_without_neighbor = view(nlist[nlist.!=nID], :)
+  balltree = BallTree(transpose(coordinates[nlist_without_neighbor, :]))
+  return nlist_without_neighbor[inrange(balltree, coordinates[nID, :], bond_horizon, true)]
 
 end
 
@@ -181,8 +181,7 @@ function compute_bond_associated_weighted_volume(nodes::Union{SubArray,Vector{In
     neighborhood_volume = sum(volume[nlist[iID]] .* bond_damage[iID] .* omega[iID])
 
     for (jID, nID) in enumerate(nlist[iID])
-      nlist_without_neighbor = view(nlist[iID][nlist[iID].!=nID], :)
-      local_nlist = find_local_neighbors(coordinates[nID, :], coordinates[nlist_without_neighbor, :], nlist_without_neighbor, bond_horizon)
+      local_nlist = find_local_neighbors(nID, coordinates, nlist[iID], bond_horizon)
       sub_bond_list = [findfirst(x -> x == elem, nlist[iID]) for elem in local_nlist]
       weighted_volume[iID][jID] = sum(volume[local_nlist] .* bond_damage[iID][sub_bond_list] .* omega[iID][sub_bond_list] / neighborhood_volume)
     end

src/Physics/Material/Material_Models/Correspondence.jl

@@ -65,8 +65,7 @@ function init_material_model(datamanager::Module, nodes::Union{SubArray,Vector{I
 
   end
   if haskey(material_parameter, "Bond Associated") && material_parameter["Bond Associated"]
-    datamanager = Bond_Associated_Correspondence.init_material_model(datamanager, nodes, material_parameter)
-    return datamanager
+    return Bond_Associated_Correspondence.init_material_model(datamanager, nodes, material_parameter)
   end
   material_parameter["Bond Associated"] = false
   return datamanager

src/Physics/Pre_calculation/Pre_Calculation_Factory.jl

@@ -69,20 +69,24 @@ function init_pre_calculation(datamanager::Module, options::Dict)
         datamanager.create_bond_field("Deformed Bond Geometry", Float64, dof)
         datamanager.create_bond_field("Deformed Bond Length", Float64, 1)
     end
+
+    if options["Deformation Gradient"]
+        datamanager.create_constant_node_field("Deformation Gradient", Float64, "Matrix", dof)
+        options["Shape Tensor"] = true
+    end
+
+    if options["Bond Associated Deformation Gradient"]
+        datamanager.create_constant_bond_field("Bond Associated Deformation Gradient", Float64, "Matrix", dof)
+        options["Bond Associated Shape Tensor"] = true
+    end
     if options["Shape Tensor"]
         datamanager.create_constant_node_field("Shape Tensor", Float64, "Matrix", dof)
         datamanager.create_constant_node_field("Inverse Shape Tensor", Float64, "Matrix", dof)
     end
-    if options["Deformation Gradient"]
-        datamanager.create_constant_node_field("Deformation Gradient", Float64, "Matrix", dof)
-    end
     if options["Bond Associated Shape Tensor"]
         datamanager.create_constant_bond_field("Bond Associated Shape Tensor", Float64, "Matrix", dof)
         datamanager.create_constant_bond_field("Inverse Bond Associated Shape Tensor", Float64, "Matrix", dof)
     end
-    if options["Bond Associated Deformation Gradient"]
-        datamanager.create_constant_bond_field("Bond Associated Deformation Gradient", Float64, "Matrix", dof)
-    end
     return datamanager
 end
 

src/Physics/Pre_calculation/bond_deformation.jl

@@ -20,12 +20,11 @@ Compute the bond deformation.
 - `datamanager`: Datamanager.
 """
 function compute(datamanager::Module, nodes::Union{SubArray,Vector{Int64}}, time::Float64)
-    dof = datamanager.get_dof()
     nlist = datamanager.get_nlist()
     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, dof, nlist, deformed_coor, deformed_bond, deformed_bond_length)
+    deformed_bond, deformed_bond_length = Geometry.bond_geometry(nodes, nlist, deformed_coor, deformed_bond, deformed_bond_length)
     return datamanager
 end
 

src/Physics/Pre_calculation/bond_deformation_gradient.jl

@@ -46,7 +46,7 @@ function compute(datamanager::Module, nodes::Union{SubArray,Vector{Int64}})
             bond_horizon = horizon[iID]
         end
         for (jID, nID) in enumerate(nlist[iID])
-            neighbor_nlist = find_local_neighbors(neighbor_coordinate[nID], coordinates[nlist[nlist.!=nID], :], nlist[iID], bond_horizon)
+            neighbor_nlist = find_local_neighbors(nID, coordinates, nlist[iID], bond_horizon)
             indices = vcat(1:jID-1, jID+1:length(nlist[iID]))
 
             deformation_gradient[iID][jID, :, :] = Geometry.bond_associated_deformation_gradient(dof, volume[neighbor_nlist], omega[neighbor_nlist], bond_damage[iID][indices], undeformed_bond[iID][indices], deformed_bond[iID][indices], bond_deformation_gradient[iID][jID, :, :])

src/Physics/Pre_calculation/bond_shape_tensor.jl

@@ -4,6 +4,7 @@
 
 module Bond_Shape_Tensor
 include("../Material/Material_Models/Bond_Associated_Correspondence.jl")
+
 using .Bond_Associated_Correspondence: find_local_neighbors
 
 export compute
@@ -26,7 +27,7 @@ function compute(datamanager::Module, nodes::Union{SubArray,Vector{Int64}})
     volume = datamanager.get_field("Volume")
     omega = datamanager.get_field("Influence Function")
     bond_damage = datamanager.get_bond_damage("NP1")
-    undeformed_bond = datamanager.get_field("Bond Geometry")
+    coordinates = datamanager.get_field("Coordinates")
     bond_shape_tensor = datamanager.get_field("Bond Associated Shape Tensor")
     inverse_bond_shape_tensor = datamanager.get_field("Inverse Bond Associated Shape Tensor")
 
@@ -42,10 +43,17 @@ function compute(datamanager::Module, nodes::Union{SubArray,Vector{Int64}})
         if isnothing(bond_horizon)
             bond_horizon = horizon[iID]
         end
+        #bond damage
         for (jID, nID) in enumerate(nlist[iID])
-            neighbor_nlist = find_local_neighbors(neighbor_coordinate[nID], coordinates[nlist[nlist.!=nID], :], nlist[iID], bond_horizon)
+
+            neighbor_nlist = find_local_neighbors(nID, coordinates, nlist[iID], bond_horizon)
+            undeformed_bond, distances = Geometry.calculate_bond_length(coor[nID, :], neighbor_nlist)
+            # TODO Bond damage is not correct and must be adapted
+            # indices are not needed for that
             indices = vcat(1:jID-1, jID+1:length(nlist[iID]))
-            shape_tensor[iID][jID, :, :], inverse_shape_tensor[iID][jID, :, :] = Geometry.bond_associated_shape_tensor(dof, volume[neighbor_nlist], omega[neighbor_nlist], bond_damage[iID][indices], undeformed_bond[iID][indices], bond_shape_tensor[iID][jID, :, :], inverse_bond_shape_tensor[iID][jID, :, :])
+
+
+            shape_tensor[iID][jID, :, :], inverse_shape_tensor[iID][jID, :, :] = Geometry.bond_associated_shape_tensor(dof, volume[neighbor_nlist], omega[neighbor_nlist], bond_damage[iID][indices], undeformed_bond, bond_shape_tensor[iID][jID, :, :], inverse_bond_shape_tensor[iID][jID, :, :])
         end
 
     end

src/Support/Parameters/parameter_handling.jl

@@ -116,6 +116,7 @@ global expected_structure = Dict(
                             "Any" => [Dict{Any,Any}(
                                     "Material Model" => [String, true],
                                     "Symmetry" => [String, false],
+                                    "Bond Associated" => [Bool, false],
                                     "Poisson's Ratio" => [Union{Float64,Int64}, false],
                                     "Young's Modulus" => [Union{Float64,Int64}, false],
                                     "Bulk Modulus" => [Union{Float64,Int64}, false],

src/Support/geometry.jl

@@ -9,13 +9,12 @@ export bond_geometry
 export shape_tensor
 
 """
-     bond_geometry(nodes::Union{SubArray,Vector{Int64}}, dof::Int64, nlist, coor, undeformed_bond, undeformed_bond_length)
+     bond_geometry(nodes::Union{SubArray,Vector{Int64}}, nlist, coor, undeformed_bond, undeformed_bond_length)
 
 Calculate bond geometries between nodes based on their coordinates.
 
 # Arguments
  - `nodes::Union{SubArray,Vector{Int64}}`: A vector of integers representing node IDs.
- - `dof::Int64`: An integer representing the degrees of freedom.
  - `nlist`: A data structure (e.g., a list or array) representing neighboring node IDs for each node.
  - `coor`: A matrix representing the coordinates of each node.
  - `undeformed_bond`: A preallocated array or data structure to store bond geometries.
@@ -39,27 +38,27 @@ 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}}, dof::Int64, nlist, coor, undeformed_bond, undeformed_bond_length)
+function bond_geometry(nodes::Union{SubArray,Vector{Int64}}, nlist::Union{SubArray,Vector{Int64}}, coor::Union{SubArray,Matrix{Float64},Matrix{Int64}}, undeformed_bond, undeformed_bond_length)
 
     for iID in nodes
-
-        # Calculate bond vector and distance
-        bond_vectors = coor[nlist[iID], :] .- coor[iID, :]'
-        distances = sqrt.(sum(bond_vectors .^ 2, dims=2))
-
-        # Check for identical point coordinates
-        if any(distances .== 0)
+        undeformed_bond[iID], undeformed_bond_length[iID] = calculate_bond_length(iID, coor, nlist[iID])
+        if any(undeformed_bond_length[iID] .== 0)
+            println()
             @error "Identical point coordinates with no distance $iID"
             return nothing
         end
-
-        undeformed_bond[iID] .= bond_vectors
-        undeformed_bond_length[iID] .= distances
-
     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 = coor[nlist, :] .- coor[iID, :]'
+    distances = sqrt.(sum(bond_vectors .^ 2, dims=2))[:]
+
+    # Check for identical point coordinates
+    return bond_vectors, distances
+end
 """
     shape_tensor(nodes::Union{SubArray, Vector{Int64}}, dof::Int64, nlist, volume, omega, bond_damage, undeformed_bond, shape_tensor, inverse_shape_tensor)
 
@@ -127,7 +126,7 @@ end
 
 
 function bond_associated_shape_tensor(dof::Int64, volume, omega, bond_damage, undeformed_bond, shape_tensor, inverse_shape_tensor)
-
+    # bond geometries -> zwischen nachbar und seinen nachbarn
     shape_tensor[:, :] = calculate_shape_tensor(shape_tensor[:, :], dof, volume, omega, bond_damage, undeformed_bond)
     try
         inverse_shape_tensor[:, :] = inv(shape_tensor[:, :])
@@ -143,7 +142,7 @@ end
 
 
 function bond_associated_deformation_gradient(dof::Int64, volume, omega, bond_damage, undeformed_bond, deformed_bond, deformation_gradient)
-
+    # bond deformation -> zwischen nachbar und seinen nachbarn
     return calculate_deformation_gradient(deformation_gradient, dof, bond_damage, deformed_bond, undeformed_bond, volume, omega)
 
 end

test/unit_tests/Core/Solver/ut_Verlet.jl

@@ -87,7 +87,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(Vector(1:nnodes), dof, nlist, coor, undeformed_bond, undeformed_bond_length)
+undeformed_bond = PeriLab.IO.Geometry.bond_geometry(Vector(1:nnodes), nlist, coor, undeformed_bond, undeformed_bond_length)
 
 blocks = [1, 1, 2, 2, 1]
 blocks = test_Data_manager.set_block_list(blocks)

test/unit_tests/IO/ut_mesh_data.jl

@@ -482,7 +482,7 @@ end
     coor = test_Data_manager.get_field("Coordinates")
     undeformed_bond = 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(Vector(1:nnodes), dof, nlist, coor, undeformed_bond, undeformed_bond_length)
+    undeformed_bond, undeformed_bond_length = PeriLab.IO.Geometry.bond_geometry(Vector(1:nnodes), nlist, coor, undeformed_bond, undeformed_bond_length)
 
     @test undeformed_bond[1][1, 1] == 1
     @test undeformed_bond[1][1, 2] == 0

test/unit_tests/Physics/Material/Material_Models/ut_Bond_Associated_Correspondence.jl

@@ -104,10 +104,10 @@ end
 
     bond_horizon::Float64 = 1
 
-    @test Bond_Associated_Correspondence.find_local_neighbors(coordinates[5, :], coordinates[nlist[nlist.!=5], :], nlist[nlist.!=5], bond_horizon) == []
+    @test Bond_Associated_Correspondence.find_local_neighbors(5, coordinates, nlist, bond_horizon) == []
     bond_horizon = 2.6
-    @test Bond_Associated_Correspondence.find_local_neighbors(coordinates[5, :], coordinates[nlist[nlist.!=5], :], nlist[nlist.!=5], bond_horizon) == [2, 3, 4]
-    @test Bond_Associated_Correspondence.find_local_neighbors(coordinates[2, :], coordinates[nlist[nlist.!=2], :], nlist[nlist.!=2], bond_horizon) == [3, 4, 5]
+    @test Bond_Associated_Correspondence.find_local_neighbors(5, coordinates, nlist, bond_horizon) == [2, 3, 4]
+    @test Bond_Associated_Correspondence.find_local_neighbors(2, coordinates, nlist, bond_horizon) == [3, 4, 5]
 end
 
 @testset "ut_compute_bond_associated_weighted_volume" begin

test/unit_tests/Physics/ut_Physics_Factory.jl

@@ -3,8 +3,8 @@
 # SPDX-License-Identifier: BSD-3-Clause
 
 include("../../../src/Physics/Physics_Factory.jl")
-# include("../../../src/Core/data_manager.jl")
-# include("../../../src/Support/Parameters/parameter_handling.jl")
+#include("../../../src/PeriLab.jl")
+#using .PeriLab
 using Test
 import .Physics
 @testset "ut_get_block_model_definition" begin
@@ -80,8 +80,8 @@ end
     test_Data_manager = Physics.init_pre_calculation(test_Data_manager, options)
     @test "Deformed Bond GeometryN" in test_Data_manager.get_all_field_keys()
     @test "Deformed Bond GeometryNP1" in test_Data_manager.get_all_field_keys()
-    @test !("Shape Tensor" in test_Data_manager.get_all_field_keys())
-    @test !("Inverse Shape Tensor" in test_Data_manager.get_all_field_keys())
+    @test "Shape Tensor" in test_Data_manager.get_all_field_keys()
+    @test "Inverse Shape Tensor" in test_Data_manager.get_all_field_keys()
     @test "Deformation Gradient" in test_Data_manager.get_all_field_keys()
     @test !("Bond Associated Shape Tensor" in test_Data_manager.get_all_field_keys())
     @test !("Bond Associated Deformation Gradient" in test_Data_manager.get_all_field_keys())
@@ -97,7 +97,7 @@ end
     @test "Bond Associated Shape Tensor" in test_Data_manager.get_all_field_keys()
     @test !("Bond Associated Deformation Gradient" in test_Data_manager.get_all_field_keys())
 
-    options = Dict("Deformed Bond Geometry" => true, "Shape Tensor" => true, "Deformation Gradient" => true, "Bond Associated Shape Tensor" => true, "Bond Associated Deformation Gradient" => true)
+    options = Dict("Deformed Bond Geometry" => true, "Shape Tensor" => false, "Deformation Gradient" => true, "Bond Associated Shape Tensor" => false, "Bond Associated Deformation Gradient" => true)
 
     test_Data_manager = Physics.init_pre_calculation(test_Data_manager, options)
     @test "Deformed Bond GeometryN" in test_Data_manager.get_all_field_keys()
@@ -106,4 +106,5 @@ end
     @test "Deformation Gradient" in test_Data_manager.get_all_field_keys()
     @test "Bond Associated Shape Tensor" in test_Data_manager.get_all_field_keys()
     @test "Bond Associated Deformation Gradient" in test_Data_manager.get_all_field_keys()
+
 end