Skip to content

Commit 21ac7f1

Browse files
CWillbergCWillberg
committed
zero energy fix
1 parent 209ffff commit 21ac7f1

25 files changed

+3084
-116
lines changed

‎src/Core/Solver/Matrix_Verlet.jl‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -255,12 +255,12 @@ function run_solver(solver_options::Dict{Any,Any},
255255
end
256256
end
257257
@timeit "compute Velocity" begin
258-
@. @views vNP1[active_nodes, :] = (1 - numerical_damping) .*
259-
vN[active_nodes, :] .+
260-
0.5 * dt .* a[active_nodes, :]
258+
@. @views vNP1[active_nodes, :] = (1 - numerical_damping) *
259+
vN[active_nodes, :] +
260+
0.5 * dt * a[active_nodes, :]
261261

262-
@. @views uNP1[active_nodes, :] = uN[active_nodes, :] .+
263-
dt .* vNP1[active_nodes, :]
262+
@. @views uNP1[active_nodes, :] = uN[active_nodes, :] +
263+
dt * vNP1[active_nodes, :]
264264
end
265265

266266
@timeit "apply BC" apply_bc_dirichlet(["Displacements", "Temperature"],
@@ -283,7 +283,7 @@ function run_solver(solver_options::Dict{Any,Any},
283283
@timeit "Force computations" begin
284284
# check if valid if volume is different
285285

286-
@. @views fNP1 = force_densities_NP1[active_nodes, :]
286+
@views fNP1 = force_densities_NP1[active_nodes, :]
287287

288288
#-= f_int(K,
289289
# vec(uNP1[active_nodes,
@@ -292,11 +292,11 @@ function run_solver(solver_options::Dict{Any,Any},
292292
vec(uNP1[active_nodes,
293293
:]), sa)
294294
@. @views fNP1 .+= external_force_densities[active_nodes,
295-
:] .+
295+
:] +
296296
external_forces[active_nodes,
297-
:] ./ volume[active_nodes]
297+
:] / volume[active_nodes]
298298

299-
@. @views forces[active_nodes, :] .= force_densities_NP1[active_nodes, :] .*
299+
@. @views forces[active_nodes, :] .= force_densities_NP1[active_nodes, :] *
300300
volume[active_nodes] +
301301
external_forces[active_nodes, :]
302302
end

‎src/Core/Solver/Matrix_linear_static.jl‎

Lines changed: 29 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -134,6 +134,10 @@ function init_solver(solver_options::Dict{Any,Any},
134134
deformed_coorN .= coor
135135
deformed_coorNP1 .= coor
136136

137+
K = Data_Manager.get_stiffness_matrix()
138+
perm = create_permutation(Data_Manager.get_nnodes(), Data_Manager.get_dof())
139+
Data_Manager.set_stiffness_matrix(K[perm, perm])
140+
137141
### for coupled thermal analysis
138142

139143
#critical_time_step::Float64 = 1.0e50
@@ -263,7 +267,7 @@ function run_solver(solver_options::Dict{Any,Any},
263267
time,
264268
step_time)
265269

266-
external_force_densities .= external_forces ./ volume # it must be delta external forces
270+
@. external_force_densities = external_forces / volume # it must be delta external forces
267271

268272
@timeit "compute_models" compute_stiff_matrix_compatible_models(block_nodes,
269273
dt,
@@ -282,14 +286,14 @@ function run_solver(solver_options::Dict{Any,Any},
282286
K = Data_Manager.get_stiffness_matrix()
283287

284288
#@views external_force_densities[active_nodes, :] += force_densities_NP1[active_nodes, :]
285-
perm = create_permutation(active_nodes, Data_Manager.get_dof()) # only active node dofs are there
286-
287-
filtered_perm = filter_and_map_bc(non_BCs, active_nodes, dof,
288-
Data_Manager.get_nnodes())
289-
290-
# @info (filtered_perm)
291-
@timeit "compute_displacements" @views compute_displacements!(K[perm, perm],
292-
filtered_perm,
289+
active_dofs = get_active_dof(active_nodes, Data_Manager.get_dof(),
290+
Data_Manager.get_nnodes())
291+
active_non_BCS = filter_and_map_bc(non_BCs, active_dofs)
292+
293+
# - force_densities_NP1 because the internal forces come from thermal analysis. For static cases they are zero
294+
@timeit "compute_displacements" @views compute_displacements!(K[active_dofs,
295+
active_dofs],
296+
active_non_BCS,
293297
uNP1[active_nodes,
294298
:],
295299
-force_densities_NP1[active_nodes,
@@ -307,11 +311,11 @@ function run_solver(solver_options::Dict{Any,Any},
307311
dt)
308312

309313
for iID in nodes
310-
@views forces[iID, :] .= force_densities_NP1[iID, :] .* volume[iID]
314+
@. @views forces[iID, :] = force_densities_NP1[iID, :] * volume[iID]
311315
end
312316

313-
velocities .= (uNP1 - uN) ./ dt
314-
@views deformed_coorNP1 .= coor .+ uNP1
317+
@. velocities = (uNP1 - uN) / dt
318+
@. @views deformed_coorNP1 = coor + uNP1
315319

316320
# @timeit "download_from_cores" Data_Manager.synch_manager(synchronise_field,
317321
# "download_from_cores")
@@ -343,28 +347,32 @@ function run_solver(solver_options::Dict{Any,Any},
343347
barrier(comm)
344348
end
345349
end
346-
Data_Manager.set_current_time(time-dt)
350+
Data_Manager.set_current_time(time - dt)
347351
return result_files
348352
end
349353

350-
function filter_and_map_bc(non_BCs, active_nodes, dof, nnodes::Int64)
351-
filtered = []
352-
for i in eachindex(active_nodes)
354+
function get_active_dof(active_nodes, dof, nnodes)
355+
active_dofs = []
356+
for iD in active_nodes
353357
for j in 1:dof
354-
if active_nodes[i] + (j - 1) * nnodes in non_BCs
355-
push!(filtered, i + (j - 1) * length(active_nodes))
356-
end
358+
push!(active_dofs, j + (iD - 1) * dof)
357359
end
358360
end
361+
return active_dofs
362+
end
359363

360-
return filtered
364+
function filter_and_map_bc(non_BCs, active_dof)
365+
return intersect(non_BCs, active_dof)
361366
end
362367

363368
function compute_matrix(nodes::AbstractVector{Int64})
364369
if length(nodes) == 0
365370
return Data_Manager.get_stiffness_matrix()
366371
end
367372
Correspondence_matrix_based.compute_model(nodes)
373+
K = Data_Manager.get_stiffness_matrix()
374+
perm = create_permutation(Data_Manager.get_nnodes(), Data_Manager.get_dof())
375+
Data_Manager.set_stiffness_matrix(K[perm, perm])
368376
end
369377
"""
370378
compute_displacements!(K, non_BCs, u, F, F_temp, K_reduced, lu_fact, temp)
@@ -403,7 +411,7 @@ function compute_displacements!(K::AbstractMatrix{Float64},
403411
# F_from_BCs = zeros(length(non_BCs))
404412
end
405413
# 3. Modified force: F_total - K_fb * u_b
406-
F_modified = F_ext[non_BCs] .- F_int[non_BCs]
414+
@views F_modified = F_ext[non_BCs] .- F_int[non_BCs]
407415
#@info non_BCs
408416
# 4. Solve for free DOFs
409417
@views vec(u)[non_BCs] .= K[non_BCs, non_BCs] \ F_modified

‎src/Core/Solver/Verlet_solver.jl‎

Lines changed: 11 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -271,8 +271,8 @@ function run_solver(solver_options::Dict{Any,Any},
271271
apply_bc_dirichlet(["Velocity"], bcs, time,
272272
step_time)
273273
@. @views uNP1[active_nodes,
274-
:] = uN[active_nodes, :] .+
275-
dt .* vNP1[active_nodes, :]
274+
:] = uN[active_nodes, :] +
275+
dt * vNP1[active_nodes, :]
276276
end
277277

278278
compute_parabolic_problems_before_model_evaluation(active_nodes,
@@ -289,7 +289,7 @@ function run_solver(solver_options::Dict{Any,Any},
289289
# all points to guarantee that the neighbors have coor as coordinates if they are not active
290290
if "Material" in solver_options["Models"]
291291
@. @views deformed_coorNP1[active_nodes,
292-
:] = coor[active_nodes, :] .+
292+
:] = coor[active_nodes, :] +
293293
uNP1[active_nodes, :]
294294
else
295295
deformed_coorNP1 = Data_Manager.get_field("Deformed Coordinates", "NP1")
@@ -337,12 +337,12 @@ function run_solver(solver_options::Dict{Any,Any},
337337
@. forces[active_nodes, :] += external_forces[active_nodes, :]
338338
@. force_densities[active_nodes,
339339
:] += external_force_densities[active_nodes,
340-
:] .+
340+
:] +
341341
external_forces[active_nodes,
342-
:] ./
342+
:] /
343343
volume[active_nodes]
344344
@. a[active_nodes,
345-
:] = force_densities[active_nodes, :] ./
345+
:] = force_densities[active_nodes, :] /
346346
lumped_mass[active_nodes] # element wise
347347

348348
active_nodes = Data_Manager.get_field("Active Nodes")
@@ -353,17 +353,14 @@ function run_solver(solver_options::Dict{Any,Any},
353353
end
354354

355355
@. @views forces[active_nodes, :] = external_forces[active_nodes, :]
356-
@. @views force_densities[active_nodes,
357-
:] += external_force_densities[active_nodes,
358-
:] .+
356+
@. @views force_densities[active_nodes, :] += external_force_densities[active_nodes,
357+
:] +
359358
external_forces[active_nodes,
360-
:] ./
359+
:] /
361360
volume[active_nodes]
362-
@. @views a[active_nodes,
363-
:] = force_densities[active_nodes, :] ./
361+
@. @views a[active_nodes, :] = force_densities[active_nodes, :] /
364362
density[active_nodes] # element wise
365-
@. @views forces[active_nodes,
366-
:] += force_densities[active_nodes, :] .*
363+
@. @views forces[active_nodes, :] += force_densities[active_nodes, :] *
367364
volume[active_nodes]
368365
end
369366

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

Lines changed: 55 additions & 72 deletions
Original file line numberDiff line numberDiff line change
@@ -408,9 +408,6 @@ end
408408
V[start_idx + (j - 1) * dof + i - 1] += scale * block[i, j]
409409
end
410410
end
411-
412-
# [Include all other functions from working version]
413-
414411
function add_zero_energy_stiff!(K::SparseMatrixCSC{Float64,Int64},
415412
active_nodes::AbstractVector{Int64},
416413
dof::Int64,
@@ -424,112 +421,98 @@ function add_zero_energy_stiff!(K::SparseMatrixCSC{Float64,Int64},
424421
I_indices = Int[]
425422
J_indices = Int[]
426423
values = Float64[]
427-
# TODO: not very memory efficient
428-
S_matrices = Vector{Matrix{Float64}}(undef, maximum(active_nodes))
429-
D_inv_X = Vector{Vector{Vector{Float64}}}(undef, maximum(active_nodes))
430424

431-
for i in active_nodes
425+
# Pre-compute α_ipj = ω_ip * V_p * (X_ip^T * D_i^{-1} * X_ij)
426+
alpha_matrices = Vector{Matrix{Float64}}(undef, maximum(active_nodes))
427+
428+
# 1. Compute all α coefficients
429+
@inbounds for i in active_nodes
432430
neighbors = nlist[i]
433431
n_neighbors = length(neighbors)
434-
@views D_inv_i = inverse_shape_tensor[i, :, :]
435-
#
436-
D_inv_X[i] = [D_inv_i * bond_geometry[i][idx] for idx in 1:n_neighbors]
432+
D_inv_i = @view inverse_shape_tensor[i, :, :]
437433

438-
S_matrices[i] = zeros(n_neighbors, n_neighbors)
439-
for p in 1:n_neighbors
440-
if bond_damage[i][p] == 0
434+
# Pre-compute D^{-1} * X_iq for all bonds
435+
D_inv_X = [D_inv_i * bond_geometry[i][q_idx] for q_idx in 1:n_neighbors]
436+
437+
alpha_matrices[i] = zeros(n_neighbors, n_neighbors)
438+
439+
for p_idx in 1:n_neighbors
440+
if bond_damage[i][p_idx] == 0.0
441441
continue
442442
end
443-
for q in 1:n_neighbors
444-
if bond_damage[i][q] == 0
445-
continue
446-
end
447-
#TODO bond geometry must be zero for some reason??
448-
S_matrices[i][p, q] = 0.0 .* dot(bond_geometry[i][p], D_inv_X[i][q])
443+
444+
p = neighbors[p_idx]
445+
V_p = volume[p]
446+
ω_ip = omega[i][p_idx] * bond_damage[i][p_idx]
447+
448+
for q_idx in 1:n_neighbors
449+
# α_ipj = ω_ip * V_p * (X_ip^T * D_i^{-1} * X_ij)
450+
alpha_matrices[i][p_idx, q_idx] = ω_ip * V_p *
451+
dot(bond_geometry[i][p_idx],
452+
D_inv_X[q_idx])
449453
end
450454
end
451455
end
452456

453-
for i in active_nodes
457+
# 2. Assemble stiffness matrix
458+
@inbounds for i in active_nodes
454459
neighbors = nlist[i]
455-
@views Z_i = zStiff[i, :, :]
456-
V_i = volume[i]
460+
Z_i = @view zStiff[i, :, :]
457461

458462
for (idx_j, j) in enumerate(neighbors)
459-
scalar_sum_j = 0.0
460-
for idx_k in eachindex(neighbors)
461-
k = neighbors[idx_k]
462-
scalar_sum_j += omega[i][idx_k] * bond_damage[i][idx_k] * volume[k] *
463-
S_matrices[i][idx_k, idx_j]
463+
if bond_damage[i][idx_j] == 0.0
464+
continue
464465
end
465466

466-
K_ji = -V_i * (1.0 - scalar_sum_j) * Z_i
467-
add_block_to_coo!(I_indices, J_indices, values, K_ji, j, i, dof)
467+
ω_ij = omega[i][idx_j] * bond_damage[i][idx_j]
468+
469+
# K^S_ij = ω_ij * Z_i * (1 - α_ijj) - Σ_{p≠j} ω_ip * Z_i * α_ijp
470+
471+
# First term: ω_ij * Z_i * (1 - α_ijj)
472+
α_ijj = alpha_matrices[i][idx_j, idx_j]
473+
K_ij = ω_ij * (1.0 - α_ijj) * Z_i
468474

469-
for (idx_k, k) in enumerate(neighbors)
470-
if k == j
475+
# Second term: subtract contributions from other neighbors
476+
for (idx_p, p) in enumerate(neighbors)
477+
if idx_p == idx_j || bond_damage[i][idx_p] == 0.0
471478
continue
472479
end
473480

474-
s_value = omega[i][idx_k] * bond_damage[i][idx_k] * volume[k] *
475-
S_matrices[i][idx_k, idx_j]
476-
K_jk = -V_i * s_value * Z_i
477-
add_block_to_coo!(I_indices, J_indices, values, K_jk, j, k, dof)
478-
end
479-
end
480-
end
481+
ω_ip = omega[i][idx_p] * bond_damage[i][idx_p]
482+
α_ijp = alpha_matrices[i][idx_j, idx_p]
481483

482-
for i in active_nodes
483-
neighbors = nlist[i]
484-
@views Z_i = zStiff[i, :, :]
485-
V_i = volume[i]
486-
487-
for (idx_j, j) in enumerate(neighbors)
488-
K_ij_total = zeros(dof, dof)
489-
490-
for (idx_bond, bond_node) in enumerate(neighbors)
491-
if idx_bond == idx_j
492-
s_jj = S_matrices[i][idx_j, idx_j]
493-
K_ij_total .-= V_i *
494-
(1.0 -
495-
omega[i][idx_j] * bond_damage[i][idx_j] * volume[j] *
496-
s_jj) * Z_i
497-
else
498-
s_kj = S_matrices[i][idx_bond, idx_j]
499-
K_ij_total .+= V_i * omega[i][idx_j] * bond_damage[i][idx_j] *
500-
volume[j] * s_kj * Z_i
501-
end
484+
K_ij .-= ω_ip * α_ijp * Z_i
502485
end
503486

504-
add_block_to_coo!(I_indices, J_indices, values, K_ij_total, i, j, dof)
487+
add_block_to_coo!(I_indices, J_indices, values, K_ij, i, j, dof)
505488
end
506489
end
507490

491+
# 3. Build sparse matrix from COO format
508492
K_stab = sparse(I_indices, J_indices, values, size(K)...)
509493

510-
for i in active_nodes
494+
# 4. Compute diagonal terms to ensure equilibrium
495+
@inbounds for i in active_nodes
511496
K_ii_block = zeros(dof, dof)
512497

513498
for d1 in 1:dof
514499
row = (i - 1) * dof + d1
515-
for j in active_nodes
516-
if i != j
517-
for d2 in 1:dof
518-
col = (j - 1) * dof + d2
519-
K_ii_block[d1, d2] -= K_stab[row, col]
520-
end
500+
for j in 1:size(K_stab, 2)
501+
if div(j - 1, dof) + 1 != i # if column node != i
502+
K_ii_block[d1, (j - 1) % dof + 1] -= K_stab[row, j]
521503
end
522504
end
523505
end
524506

525-
for d1 in 1:dof
526-
for d2 in 1:dof
527-
row = (i - 1) * dof + d1
528-
col = (i - 1) * dof + d2
529-
K[row, col] += K_ii_block[d1, d2]
530-
end
507+
# Add diagonal block to K
508+
for d1 in 1:dof, d2 in 1:dof
509+
row = (i - 1) * dof + d1
510+
col = (i - 1) * dof + d2
511+
K[row, col] += K_ii_block[d1, d2]
531512
end
532513
end
514+
515+
# 5. Add off-diagonal stabilization terms
533516
K .+= K_stab
534517

535518
return nothing

0 commit comments

Comments
 (0)