some fixes

Author: CWillberg

src/Core/Solver/Matrix_linear_static.jl

@@ -122,7 +122,7 @@ function init_solver(solver_options::Dict{Any,Any},
 
     for (block, nodes) in pairs(block_nodes)
         model_param = Data_Manager.get_properties(block, "Material Model")
-        Correspondence_matrix_based.init_model(nodes, model_param)
+        Correspondence_matrix_based.init_model(nodes, model_param, block)
     end
     @timeit "init_matrix" Correspondence_matrix_based.init_matrix()
 
@@ -217,7 +217,7 @@ function run_solver(solver_options::Dict{Any,Any},
                                      nodes)
 
     if matrix_update
-        @timeit "update stiffness matrix" K = compute_matrix(active_nodes)
+        @timeit "update stiffness matrix" K=compute_matrix(active_nodes)
     else
         K = Data_Manager.get_stiffness_matrix()
     end
@@ -227,13 +227,16 @@ function run_solver(solver_options::Dict{Any,Any},
 
             # reshape
             # All 2 time steps fields must be in the loop, because switch NP1 to N changes the adresses
-            uN = Data_Manager.get_field("Displacements", "N")
-            uNP1 = Data_Manager.get_field("Displacements", "NP1")
-            velocities = Data_Manager.get_field("Velocity", "NP1")
-            deformed_coorNP1 = Data_Manager.get_field("Deformed Coordinates", "NP1")
-            forces = Data_Manager.get_field("Forces", "NP1")
-            force_densities_N = Data_Manager.get_field("Force Densities", "N")
-            force_densities_NP1 = Data_Manager.get_field("Force Densities", "NP1")
+            uN::NodeVectorField{Float64} = Data_Manager.get_field("Displacements", "N")
+            uNP1::NodeVectorField{Float64} = Data_Manager.get_field("Displacements", "NP1")
+            velocities::NodeVectorField{Float64} = Data_Manager.get_field("Velocity", "NP1")
+            deformed_coorNP1::NodeVectorField{Float64} = Data_Manager.get_field("Deformed Coordinates",
+                                                                                "NP1")
+            forces::NodeVectorField{Float64} = Data_Manager.get_field("Forces", "NP1")
+            force_densities_N::NodeVectorField{Float64} = Data_Manager.get_field("Force Densities",
+                                                                                 "N")
+            force_densities_NP1::NodeVectorField{Float64} = Data_Manager.get_field("Force Densities",
+                                                                                   "NP1")
 
             non_BCs = Data_Manager.get_bc_free_dof()
             #force_densities_NP1[active_nodes, :] .= force_densities_N[active_nodes, :]
@@ -258,6 +261,7 @@ function run_solver(solver_options::Dict{Any,Any},
                                step_time)
 
             external_force_densities .= external_forces ./ volume # it must be delta external forces
+            deformed_coorNP1 .= coor .+ uNP1
 
             @timeit "compute_models" compute_stiff_matrix_compatible_models(block_nodes,
                                                                             dt,
@@ -272,7 +276,7 @@ function run_solver(solver_options::Dict{Any,Any},
                                              nodes)
 
             if matrix_update
-                @timeit "update stiffness matrix" K = compute_matrix(active_nodes)
+                @timeit "update stiffness matrix" K=compute_matrix(active_nodes)
 
             else
                 K = Data_Manager.get_stiffness_matrix()
@@ -346,8 +350,8 @@ function filter_and_map_bc(non_BCs, active_nodes, dof, nnodes::Int64)
     filtered = []
     for i in eachindex(active_nodes)
         for j in 1:dof
-            if active_nodes[i] + (j - 1)*nnodes in non_BCs
-                push!(filtered, i+(j - 1) * length(active_nodes))
+            if active_nodes[i] + (j - 1) * nnodes in non_BCs
+                push!(filtered, i + (j - 1) * length(active_nodes))
             end
         end
     end
@@ -356,7 +360,7 @@ function filter_and_map_bc(non_BCs, active_nodes, dof, nnodes::Int64)
 end
 
 function compute_matrix(nodes::AbstractVector{Int64})
-    if length(nodes)==0
+    if length(nodes) == 0
         return Data_Manager.get_stiffness_matrix()
     end
     Bond_Deformation.compute(nodes,

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

@@ -12,7 +12,7 @@ using StaticArrays: @MMatrix
 using ...Data_Manager
 using ...Material_Basis: get_Hooke_matrix
 using ....Helpers: get_fourth_order, progress_bar
-using ...Zero_Energy_Control: create_zero_energy_mode_stiffness!
+using ...Zero_Energy_Control
 
 export init_model
 export add_zero_energy_stiff!
@@ -129,36 +129,45 @@ function compute_linearized_operator(CB_ik::AbstractArray{Float64,3},
         factor = omega_ij * V_k * omega_ik
 
         @inbounds begin
-            K_ijk[1, 1] = factor * (CB_ik[1, 1, 1]*DX_1 + CB_ik[1, 2, 1]*DX_2)
-            K_ijk[1, 2] = factor * (CB_ik[1, 1, 2]*DX_1 + CB_ik[1, 2, 2]*DX_2)
-            K_ijk[2, 1] = factor * (CB_ik[2, 1, 1]*DX_1 + CB_ik[2, 2, 1]*DX_2)
-            K_ijk[2, 2] = factor * (CB_ik[2, 1, 2]*DX_1 + CB_ik[2, 2, 2]*DX_2)
+            K_ijk[1, 1] = factor * (CB_ik[1, 1, 1] * DX_1 + CB_ik[1, 2, 1] * DX_2)
+            K_ijk[1, 2] = factor * (CB_ik[1, 1, 2] * DX_1 + CB_ik[1, 2, 2] * DX_2)
+            K_ijk[2, 1] = factor * (CB_ik[2, 1, 1] * DX_1 + CB_ik[2, 2, 1] * DX_2)
+            K_ijk[2, 2] = factor * (CB_ik[2, 1, 2] * DX_1 + CB_ik[2, 2, 2] * DX_2)
         end
     elseif dof == 3
-        DX_1 = D_inv_i[1, 1]*X_ij[1] + D_inv_i[1, 2]*X_ij[2] + D_inv_i[1, 3]*X_ij[3]
-        DX_2 = D_inv_i[2, 1]*X_ij[1] + D_inv_i[2, 2]*X_ij[2] + D_inv_i[2, 3]*X_ij[3]
-        DX_3 = D_inv_i[3, 1]*X_ij[1] + D_inv_i[3, 2]*X_ij[2] + D_inv_i[3, 3]*X_ij[3]
+        DX_1 = D_inv_i[1, 1] * X_ij[1] + D_inv_i[1, 2] * X_ij[2] + D_inv_i[1, 3] * X_ij[3]
+        DX_2 = D_inv_i[2, 1] * X_ij[1] + D_inv_i[2, 2] * X_ij[2] + D_inv_i[2, 3] * X_ij[3]
+        DX_3 = D_inv_i[3, 1] * X_ij[1] + D_inv_i[3, 2] * X_ij[2] + D_inv_i[3, 3] * X_ij[3]
         factor = omega_ij * V_k * omega_ik
 
         @inbounds begin
             K_ijk[1, 1] = factor *
-                          (CB_ik[1, 1, 1]*DX_1 + CB_ik[1, 2, 1]*DX_2 + CB_ik[1, 3, 1]*DX_3)
+                          (CB_ik[1, 1, 1] * DX_1 + CB_ik[1, 2, 1] * DX_2 +
+                           CB_ik[1, 3, 1] * DX_3)
             K_ijk[1, 2] = factor *
-                          (CB_ik[1, 1, 2]*DX_1 + CB_ik[1, 2, 2]*DX_2 + CB_ik[1, 3, 2]*DX_3)
+                          (CB_ik[1, 1, 2] * DX_1 + CB_ik[1, 2, 2] * DX_2 +
+                           CB_ik[1, 3, 2] * DX_3)
             K_ijk[1, 3] = factor *
-                          (CB_ik[1, 1, 3]*DX_1 + CB_ik[1, 2, 3]*DX_2 + CB_ik[1, 3, 3]*DX_3)
+                          (CB_ik[1, 1, 3] * DX_1 + CB_ik[1, 2, 3] * DX_2 +
+                           CB_ik[1, 3, 3] * DX_3)
             K_ijk[2, 1] = factor *
-                          (CB_ik[2, 1, 1]*DX_1 + CB_ik[2, 2, 1]*DX_2 + CB_ik[2, 3, 1]*DX_3)
+                          (CB_ik[2, 1, 1] * DX_1 + CB_ik[2, 2, 1] * DX_2 +
+                           CB_ik[2, 3, 1] * DX_3)
             K_ijk[2, 2] = factor *
-                          (CB_ik[2, 1, 2]*DX_1 + CB_ik[2, 2, 2]*DX_2 + CB_ik[2, 3, 2]*DX_3)
+                          (CB_ik[2, 1, 2] * DX_1 + CB_ik[2, 2, 2] * DX_2 +
+                           CB_ik[2, 3, 2] * DX_3)
             K_ijk[2, 3] = factor *
-                          (CB_ik[2, 1, 3]*DX_1 + CB_ik[2, 2, 3]*DX_2 + CB_ik[2, 3, 3]*DX_3)
+                          (CB_ik[2, 1, 3] * DX_1 + CB_ik[2, 2, 3] * DX_2 +
+                           CB_ik[2, 3, 3] * DX_3)
             K_ijk[3, 1] = factor *
-                          (CB_ik[3, 1, 1]*DX_1 + CB_ik[3, 2, 1]*DX_2 + CB_ik[3, 3, 1]*DX_3)
+                          (CB_ik[3, 1, 1] * DX_1 + CB_ik[3, 2, 1] * DX_2 +
+                           CB_ik[3, 3, 1] * DX_3)
             K_ijk[3, 2] = factor *
-                          (CB_ik[3, 1, 2]*DX_1 + CB_ik[3, 2, 2]*DX_2 + CB_ik[3, 3, 2]*DX_3)
+                          (CB_ik[3, 1, 2] * DX_1 + CB_ik[3, 2, 2] * DX_2 +
+                           CB_ik[3, 3, 2] * DX_3)
             K_ijk[3, 3] = factor *
-                          (CB_ik[3, 1, 3]*DX_1 + CB_ik[3, 2, 3]*DX_2 + CB_ik[3, 3, 3]*DX_3)
+                          (CB_ik[3, 1, 3] * DX_1 + CB_ik[3, 2, 3] * DX_2 +
+                           CB_ik[3, 3, 3] * DX_3)
         end
     end
 end
@@ -224,19 +233,19 @@ function create_mapping_structure(nodes::AbstractVector{Int64},
         mapping[i][1] = idx
         for n in 1:dof
             for m in 1:dof
-                I[idx] = (i-1)*dof + m
-                J[idx] = (i-1)*dof + n
+                I[idx] = (i - 1) * dof + m
+                J[idx] = (i - 1) * dof + n
                 idx += 1
             end
         end
 
         # Blocks (i, k) for k ∈ nlist[i]
         for (k_idx, k) in enumerate(ni)
-            mapping[i][k_idx+1] = idx
+            mapping[i][k_idx + 1] = idx
             for n in 1:dof
                 for m in 1:dof
-                    I[idx] = (i-1)*dof + m
-                    J[idx] = (k-1)*dof + n
+                    I[idx] = (i - 1) * dof + m
+                    J[idx] = (k - 1) * dof + n
                     idx += 1
                 end
             end
@@ -322,7 +331,7 @@ function assemble_stiffness(nodes::AbstractVector{Int64},
             mapping_i = mapping[i]
 
             for (j_idx, j) in enumerate(ni)
-                omega_ij = omega[i][j_idx]*bond_damage[i][j_idx]
+                omega_ij = omega[i][j_idx] * bond_damage[i][j_idx]
                 V_j = volume[j]
                 X_ij = bond_geometry[i][j_idx]
 
@@ -331,7 +340,7 @@ function assemble_stiffness(nodes::AbstractVector{Int64},
                 mapping_j = mapping[j]
 
                 for (k_idx, k) in enumerate(ni)
-                    omega_ik = omega[i][k_idx]*bond_damage[i][k_idx]
+                    omega_ik = omega[i][k_idx] * bond_damage[i][k_idx]
                     V_k = volume[k]
 
                     compute_linearized_operator(@view(CB_k[k_idx, :, :, :]),
@@ -341,7 +350,7 @@ function assemble_stiffness(nodes::AbstractVector{Int64},
                     # Row j entries
                     i_pos_in_j = neighbor_pos_j[i]
                     if i_pos_in_j != 0
-                        add_block_scaled!(V_thread, mapping_j[i_pos_in_j+1],
+                        add_block_scaled!(V_thread, mapping_j[i_pos_in_j + 1],
                                           buffers.K_block, V_i, dof)
                     end
 
@@ -350,14 +359,14 @@ function assemble_stiffness(nodes::AbstractVector{Int64},
                     else
                         k_pos_in_j = neighbor_pos_j[k]
                         if k_pos_in_j != 0
-                            add_block_scaled!(V_thread, mapping_j[k_pos_in_j+1],
+                            add_block_scaled!(V_thread, mapping_j[k_pos_in_j + 1],
                                               buffers.K_block, -V_i, dof)
                         end
                     end
 
                     # Row i entries
                     add_block_scaled!(V_thread, i_diag_idx, buffers.K_block, -V_j, dof)
-                    add_block_scaled!(V_thread, mapping_i[k_idx+1], buffers.K_block, V_j,
+                    add_block_scaled!(V_thread, mapping_i[k_idx + 1], buffers.K_block, V_j,
                                       dof)
                 end
             end
@@ -396,7 +405,7 @@ end
 
 @inline function add_block_scaled!(V, start_idx, block, scale, dof)
     @turbo for i in 1:dof, j in 1:dof
-        V[start_idx+(j-1)*dof+i-1] += scale * block[i, j]
+        V[start_idx + (j - 1) * dof + i - 1] += scale * block[i, j]
     end
 end
 
@@ -500,11 +509,11 @@ function add_zero_energy_stiff!(K::SparseMatrixCSC{Float64,Int64},
         K_ii_block = zeros(dof, dof)
 
         for d1 in 1:dof
-            row = (i-1)*dof + d1
+            row = (i - 1) * dof + d1
             for j in active_nodes
                 if i != j
                     for d2 in 1:dof
-                        col = (j-1)*dof + d2
+                        col = (j - 1) * dof + d2
                         K_ii_block[d1, d2] -= K_stab[row, col]
                     end
                 end
@@ -513,8 +522,8 @@ function add_zero_energy_stiff!(K::SparseMatrixCSC{Float64,Int64},
 
         for d1 in 1:dof
             for d2 in 1:dof
-                row = (i-1)*dof + d1
-                col = (i-1)*dof + d2
+                row = (i - 1) * dof + d1
+                col = (i - 1) * dof + d2
                 K[row, col] += K_ii_block[d1, d2]
             end
         end
@@ -534,8 +543,8 @@ function add_block_to_coo!(I_indices::Vector{Int},
     for d1 in 1:dof
         for d2 in 1:dof
             if abs(block[d1, d2]) > 1e-14
-                push!(I_indices, (i-1)*dof + d1)
-                push!(J_indices, (j-1)*dof + d2)
+                push!(I_indices, (i - 1) * dof + d1)
+                push!(J_indices, (j - 1) * dof + d2)
                 push!(values, block[d1, d2])
             end
         end
@@ -548,15 +557,15 @@ function compute_bond_force(bond_force::Vector{Vector{Vector{Float64}}},
                             dof::Int64)
     @inbounds for i in nodes
         ni = nlist[i]
-        i_offset = (i-1)*dof
+        i_offset = (i - 1) * dof
 
         for (j_idx, j) in enumerate(ni)
-            j_offset = (j-1)*dof
+            j_offset = (j - 1) * dof
             f_ij = zeros(dof)
 
             for m in 1:dof
                 for n in 1:dof
-                    f_ij[m] += K[i_offset+m, j_offset+n] * (u[j, n] - u[i, n])
+                    f_ij[m] += K[i_offset + m, j_offset + n] * (u[j, n] - u[i, n])
                 end
             end
 
@@ -566,8 +575,9 @@ function compute_bond_force(bond_force::Vector{Vector{Vector{Float64}}},
 end
 
 function init_model(nodes::AbstractVector{Int64},
-                    material_parameter::Dict)
-    if Data_Manager.get_max_rank()>1
+                    material_parameter::Dict, block_id::Int64)
+    Zero_Energy_Control.init_model(nodes, material_parameter, block_id)
+    if Data_Manager.get_max_rank() > 1
         @error "Correspondence matrix based not implemented for parallel runs."
     end
 end
@@ -581,9 +591,10 @@ function init_matrix()
     bond_geometry = Data_Manager.get_field("Bond Geometry")
     inverse_shape_tensor = Data_Manager.get_field("Inverse Shape Tensor")
     nlist = Data_Manager.get_nlist()
-    number_of_neighbors=Data_Manager.get_field("Number of Neighbors")
+    number_of_neighbors = Data_Manager.get_field("Number of Neighbors")
     volume = Data_Manager.get_field("Volume")
     omega = Data_Manager.get_field("Influence Function")
+    # using the elasticity matrix allows the introduction of heterogeneous materials
     C_voigt = Data_Manager.get_field("Elasticity Matrix")
     bond_geometry_N = Data_Manager.get_field("Deformed Bond Geometry", "N")
     bond_damage = Data_Manager.get_field("Bond Damage", "NP1")
@@ -604,23 +615,26 @@ function init_matrix()
 
     Data_Manager.init_stiffness_matrix(index_x, index_y, vals, total_dof)
     K_sparse = Data_Manager.get_stiffness_matrix()
-
-    if Data_Manager.get_analysis_model("Zero Energy Control Model", 1) != "Global"
-        @warn "Global Energy Control Model is not active for block 1. Must be for all blocks, when used."
+    if haskey(Data_Manager.get_properties(1, "Material Model"), "Zero Energy Control")
+        if Data_Manager.get_properties(1, "Material Model")["Zero Energy Control"] ==
+           "Global"
+            Zero_Energy_Control.create_zero_energy_mode_stiffness!(nodes, dof, C_voigt,
+                                                                   inverse_shape_tensor,
+                                                                   zStiff)
+
+            add_zero_energy_stiff!(K_sparse,
+                                   nodes,
+                                   dof,
+                                   zStiff,
+                                   inverse_shape_tensor,
+                                   nlist,
+                                   volume,
+                                   bond_geometry_N,
+                                   bond_damage,
+                                   omega)
+        end
     else
-        create_zero_energy_mode_stiffness!(nodes, dof, C_voigt,
-                                           inverse_shape_tensor, zStiff)
-
-        add_zero_energy_stiff!(K_sparse,
-                               nodes,
-                               dof,
-                               zStiff,
-                               inverse_shape_tensor,
-                               nlist,
-                               volume,
-                               bond_geometry_N,
-                               bond_damage,
-                               omega)
+        @warn "Global Energy Control Model is not active for block 1. Must be for all blocks, when used."
     end
 end
 
@@ -653,8 +667,8 @@ function compute_model(nodes::AbstractVector{Int64})
 
     Data_Manager.init_stiffness_matrix(index_x, index_y, vals, total_dof)
 
-    create_zero_energy_mode_stiffness!(nodes, dof, C_voigt,
-                                       inverse_shape_tensor, zStiff)
+    Zero_Energy_Control.create_zero_energy_mode_stiffness!(nodes, dof, C_voigt,
+                                                           inverse_shape_tensor, zStiff)
 
     K_sparse = Data_Manager.get_stiffness_matrix()
     add_zero_energy_stiff!(K_sparse,

src/Models/Material/Material_Models/Zero_Energy_Control/Zero_Energy_Control.jl

@@ -26,8 +26,7 @@ function init_model(nodes::AbstractVector{Int64}, material_parameter::Dict, bloc
         Data_Manager.set_model_module(zero_energy_model, mod)
         mod.init_model(nodes, material_parameter)
     else
-        Data_Manager.set_analysis_model("Zero Energy Control Model", block, [])
-
+        Data_Manager.set_analysis_model("Zero Energy Control Model", block, "")
         @warn "No zero energy control activated for corresponcence in block $block. This might cause errors."
     end
 end
@@ -57,8 +56,7 @@ function create_zero_energy_mode_stiffness!(nodes::AbstractVector{Int64},
                                             Kinv::Array{Float64,3},
                                             zStiff::Array{Float64,3})
     zero_energy_model = Data_Manager.get_analysis_model("Zero Energy Control Model", 1)
-    mod = Data_Manager.get_model_module(zero_energy_model)
-    return mod.create_zero_energy_mode_stiffness!(nodes, dof, CVoigt,
-                                                  Kinv, zStiff)
+    mod = Data_Manager.get_model_module(zero_energy_model[1])
+    return mod.create_zero_energy_mode_stiffness!(nodes, dof, CVoigt, Kinv, zStiff)
 end
 end

test/fullscale_tests/test_linear_static_matrix_based_solver/linear_static_matrix_based_solver_2D.yaml

@@ -22,11 +22,13 @@ PeriLab:
         Symmetry: "isotropic plane stress"
         Young's Modulus: 209695.43147208123
         Poisson's Ratio: 0.2944
+        Zero Energy Control: "Global"
       Mat_2:
         Material Model: "Correspondence Elastic"
         Symmetry: "isotropic plane stress"
         Young's Modulus: 209695.43147208123
         Poisson's Ratio: 0.2944
+        Zero Energy Control: "Global"
   Blocks:
     block_1:
       Block ID: 1