Skip to content

Commit

Permalink
Rework traction computation for clarity
Browse files Browse the repository at this point in the history
  • Loading branch information
@lxmota
lxmota committed Dec 7, 2024
1 parent 6110424 commit 3880590
Show file tree
Hide file tree
Showing 2 changed files with 30 additions and 35 deletions.
55 changes: 23 additions & 32 deletions src/ics_bcs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -239,6 +239,7 @@ function apply_sm_schwarz_coupling_dirichlet(model::SolidMechanics, bc::Coupling
point_posn = elem_posn * N
point_velo = elem_velo * N
point_acce = elem_acce * N
@debug "Applying Schwarz DBC as $point_posn"
model.current[:, node_index] = point_posn
model.velocity[:, node_index] = point_velo
model.acceleration[:, node_index] = point_acce
Expand All @@ -253,7 +254,8 @@ function apply_sm_schwarz_coupling_neumann(model::SolidMechanics, bc::CouplingSc
num_local_nodes = length(local_to_global_map)
for local_node ∈ 1:num_local_nodes
global_node = local_to_global_map[local_node]
node_tractions = schwarz_tractions[3*local_node-2:3*local_node]
node_tractions = schwarz_tractions[:, local_node]
@debug "Applying Schwarz NBC as $node_tractions"
model.boundary_force[3*global_node-2:3*global_node] += node_tractions
end
end
Expand Down Expand Up @@ -400,7 +402,7 @@ function apply_sm_schwarz_contact_neumann(model::SolidMechanics, bc::SMContactSc
num_local_nodes = length(local_to_global_map)
for local_node ∈ 1:num_local_nodes
global_node = local_to_global_map[local_node]
node_tractions = schwarz_tractions[3*local_node-2:3*local_node]
node_tractions = schwarz_tractions[:, local_node]
normal = normals[:, local_node]
model.boundary_force[3*global_node-2:3*global_node] += transfer_normal_component(
node_tractions,
Expand All @@ -411,18 +413,18 @@ function apply_sm_schwarz_contact_neumann(model::SolidMechanics, bc::SMContactSc
end


function reduce_traction(
function local_traction_from_global_force(
mesh::ExodusDatabase,
side_set_id::Integer,
global_traction::Vector{Float64},
global_force::Vector{Float64},
)
local_to_global_map = get_side_set_local_to_global_map(mesh, side_set_id)
num_local_nodes = length(local_to_global_map)
local_traction = zeros(3 * num_local_nodes)
local_traction = zeros(3, num_local_nodes)
for local_node ∈ 1:num_local_nodes
global_node = local_to_global_map[local_node]
local_traction[3*local_node-2:3*local_node] =
global_traction[3*global_node-2:3*global_node]
local_traction[:, local_node] =
global_force[3*global_node-2:3*global_node]
end
return local_traction
end
Expand All @@ -449,41 +451,30 @@ end
function get_dst_traction(bc::SMContactSchwarzBC)
src_mesh = bc.coupled_subsim.model.mesh
src_side_set_id = bc.coupled_side_set_id
src_global_traction = -bc.coupled_subsim.model.internal_force
src_local_traction = reduce_traction(src_mesh, src_side_set_id, src_global_traction)
src_traction_x = src_local_traction[1:3:end]
src_traction_y = src_local_traction[2:3:end]
src_traction_z = src_local_traction[3:3:end]
dst_traction_x = bc.transfer_operator * src_traction_x
dst_traction_y = bc.transfer_operator * src_traction_y
dst_traction_z = bc.transfer_operator * src_traction_z
dst_traction = zeros(3 * length(dst_traction_x))
dst_traction[1:3:end] = dst_traction_x
dst_traction[2:3:end] = dst_traction_y
dst_traction[3:3:end] = dst_traction_z
src_global_force = -bc.coupled_subsim.model.internal_force
src_local_traction = local_traction_from_global_force(src_mesh, src_side_set_id, src_global_force)
num_dst_nodes = size(bc.transfer_operator, 1)
dst_traction = zeros(3, num_dst_nodes)
dst_traction[1, :] = bc.transfer_operator * src_local_traction[1, :]
dst_traction[2, :] = bc.transfer_operator * src_local_traction[2, :]
dst_traction[3, :] = bc.transfer_operator * src_local_traction[3, :]
return dst_traction
end

function get_dst_traction(bc::SMNonOverlapSchwarzBC)
src_mesh = bc.coupled_subsim.model.mesh
src_side_set_id = bc.coupled_side_set_id
src_global_traction = -bc.coupled_subsim.model.internal_force
src_local_traction = reduce_traction(src_mesh, src_side_set_id, src_global_traction)
src_traction_x = src_local_traction[1:3:end]
src_traction_y = src_local_traction[2:3:end]
src_traction_z = src_local_traction[3:3:end]
src_global_force = -bc.coupled_subsim.model.internal_force
src_local_traction = local_traction_from_global_force(src_mesh, src_side_set_id, src_global_force)
compute_transfer_operator(bc.coupled_subsim.model, bc)
dst_traction_x = bc.transfer_operator * src_traction_x
dst_traction_y = bc.transfer_operator * src_traction_y
dst_traction_z = bc.transfer_operator * src_traction_z
dst_traction = zeros(3 * length(dst_traction_x))
dst_traction[1:3:end] = dst_traction_x
dst_traction[2:3:end] = dst_traction_y
dst_traction[3:3:end] = dst_traction_z
num_dst_nodes = size(bc.transfer_operator, 1)
dst_traction = zeros(3, num_dst_nodes)
dst_traction[1, :] = bc.transfer_operator * src_local_traction[1, :]
dst_traction[2, :] = bc.transfer_operator * src_local_traction[2, :]
dst_traction[3, :] = bc.transfer_operator * src_local_traction[3, :]
return dst_traction
end


function node_set_id_from_name(node_set_name::String, mesh::ExodusDatabase)
node_set_names = Exodus.read_names(mesh, NodeSet)
num_names = length(node_set_names)
Expand Down
10 changes: 7 additions & 3 deletions src/schwarz.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -195,8 +195,12 @@ function subcycle(sim::MultiDomainSimulation, is_schwarz::Bool)
break
end
subsim.model.time = subsim.integrator.time
@debug "Subdomain : $(subsim.name), Time : $(subsim.model.time)"
@debug "Before applying BCs"
apply_bcs(subsim)
@debug "After applying BCs and before advance"
advance(subsim)
@debug "After advance"
if subsim.failed == true
sim.failed = true
return
Expand Down Expand Up @@ -356,14 +360,14 @@ function check_compression(
)
compression_tol = 0.0
compression = false
reactions = get_dst_traction(bc)
nodal_reactions = get_dst_traction(bc)
normals = compute_normal(mesh, bc.side_set_id, model)
local_to_global_map = get_side_set_local_to_global_map(mesh, bc.side_set_id)
num_local_nodes = length(local_to_global_map)
for local_node ∈ 1:num_local_nodes
reaction_node = reactions[3*local_node-2:3*local_node]
nodal_reaction = nodal_reactions[:, local_node]
normal = normals[:, local_node]
normal_traction = dot(reaction_node, normal)
normal_traction = dot(nodal_reaction, normal)
compressive_traction = normal_traction ≤ compression_tol
if compressive_traction == true
compression = true
Expand Down

0 comments on commit 3880590

Please sign in to comment.