From 989c09480af2dba6c044274bfee0475f72ddfc81 Mon Sep 17 00:00:00 2001
From: Alejandro Mota <amota@sandia.gov>
Date: Fri, 31 Jan 2025 14:29:26 -0800
Subject: [PATCH] Run formatter
---
src/ems/cloud.jl | 8 +-
src/ics_bcs.jl | 193 ++++++++++++++++++++---------------------
src/model.jl | 16 +---
src/simulation.jl | 6 +-
src/solver.jl | 82 ++++++++---------
src/time_integrator.jl | 134 ++++++++++++++--------------
6 files changed, 212 insertions(+), 227 deletions(-)
diff --git a/src/ems/cloud.jl b/src/ems/cloud.jl
index b49a9c7..b28b60e 100644
--- a/src/ems/cloud.jl
+++ b/src/ems/cloud.jl
@@ -245,7 +245,7 @@ function compute_energy_force!(potential::SymmetricSethHill, points::Matrix{Floa
γ = r / h
κ = γ^m
α = κ - 1.0
- β = 1.0 / κ - 1.0
+ β = 1.0 / κ - 1.0
energy += (0.5 * k / m / m) * (α * α + β * β)
f_ij = (-k / m) * (α * κ / γ - β / κ / γ) * (r_ij / r)
force[:, i] .-= f_ij
@@ -382,9 +382,9 @@ function write_points_vtk(filename::String, points::Matrix{Float64})
if size(points, 1) != 3
error("The input matrix must have size 3xN, where N is the number of points.")
end
-
+
N = size(points, 2) # Number of points
-
+
# Open the file for writing
open(filename, "w") do io
# Write VTK header
@@ -393,7 +393,7 @@ function write_points_vtk(filename::String, points::Matrix{Float64})
println(io, "ASCII")
println(io, "DATASET POLYDATA")
println(io, "POINTS $N float")
-
+
# Write point data
for i in 1:N
println(io, join(points[:, i], " "))
diff --git a/src/ics_bcs.jl b/src/ics_bcs.jl
index 64abe95..cf50dc8 100644
--- a/src/ics_bcs.jl
+++ b/src/ics_bcs.jl
@@ -88,9 +88,9 @@ function SMContactSchwarzBC(
transfer_operator = Matrix{Float64}(undef, 0, 0)
rotation_matrix = I(3)
active_contact = false
- swap_bcs = false
+ swap_bcs = false
if haskey(bc_params, "swap BC types") == true
- swap_bcs = bc_params["swap BC types"]
+ swap_bcs = bc_params["swap BC types"]
end
SMContactSchwarzBC(
side_set_name,
@@ -143,10 +143,10 @@ function SMCouplingSchwarzBC(
_, _, side_set_node_indices = get_side_set_local_from_global_map(input_mesh, side_set_id)
coupled_block_name = bc_params["source block"]
if typeof(coupled_subsim.model) == LinearOpInfRom
- coupled_mesh = coupled_subsim.model.fom_model.mesh
+ coupled_mesh = coupled_subsim.model.fom_model.mesh
else
- coupled_mesh = coupled_subsim.model.mesh
- end
+ coupled_mesh = coupled_subsim.model.mesh
+ end
coupled_block_id = block_id_from_name(coupled_block_name, coupled_mesh)
element_type = Exodus.read_block_parameters(coupled_mesh, coupled_block_id)[1]
coupled_side_set_name = bc_params["source side set"]
@@ -170,7 +170,7 @@ function SMCouplingSchwarzBC(
push!(interpolation_function_values, N)
end
is_dirichlet = true
- swap_bcs = false
+ swap_bcs = false
if bc_type == "Schwarz overlap"
SMOverlapSchwarzBC(
side_set_name,
@@ -185,14 +185,14 @@ function SMCouplingSchwarzBC(
elseif bc_type == "Schwarz nonoverlap"
if haskey(bc_params, "default BC type") == true
default_bc_type = bc_params["default BC type"]
- if default_bc_type == "Dirichlet"
+ if default_bc_type == "Dirichlet"
is_dirichlet = true
- elseif default_bc_type == "Neumann"
+ elseif default_bc_type == "Neumann"
is_dirichlet = false
- else
- error("Invalid string for 'default BC type'! Valid options are 'Dirichlet' and 'Neumann'")
- end
- end
+ else
+ error("Invalid string for 'default BC type'! Valid options are 'Dirichlet' and 'Neumann'")
+ end
+ end
if haskey(bc_params, "swap BC types") == true
swap_bcs = bc_params["swap BC types"]
end
@@ -204,7 +204,7 @@ function SMCouplingSchwarzBC(
coupled_subsim,
subsim,
coupled_side_set_id,
- is_dirichlet,
+ is_dirichlet,
swap_bcs
)
else
@@ -214,30 +214,30 @@ end
function apply_bc(model::LinearOpInfRom, bc::SMDirichletBC)
model.fom_model.time = model.time
- apply_bc(model.fom_model,bc)
+ apply_bc(model.fom_model, bc)
bc_vector = zeros(0)
for node_index ∈ bc.node_set_node_indices
dof_index = 3 * (node_index - 1) + bc.offset
- disp_val = model.fom_model.current[bc.offset,node_index] - model.fom_model.reference[bc.offset, node_index]
- push!(bc_vector,disp_val)
+ disp_val = model.fom_model.current[bc.offset, node_index] - model.fom_model.reference[bc.offset, node_index]
+ push!(bc_vector, disp_val)
end
offset = bc.offset
if offset == 1
- offset_name = "x"
+ offset_name = "x"
end
if offset == 2
- offset_name = "y"
+ offset_name = "y"
end
if offset == 3
- offset_name = "z"
+ offset_name = "z"
end
- op_name = "B_"*bc.node_set_name * "-" * offset_name
- bc_operator = model.opinf_rom[op_name]
+ op_name = "B_" * bc.node_set_name * "-" * offset_name
+ bc_operator = model.opinf_rom[op_name]
# SM Dirichlet BC are only defined on a single x,y,z
- model.reduced_boundary_forcing[:] += bc_operator[1,:,:] * bc_vector
- end
+ model.reduced_boundary_forcing[:] += bc_operator[1, :, :] * bc_vector
+end
function apply_bc(model::SolidMechanics, bc::SMDirichletBC)
for node_index ∈ bc.node_set_node_indices
@@ -332,7 +332,7 @@ function apply_bc(model::SolidMechanics, bc::SMDirichletInclined)
model.velocity[:, node_index] = original_velocity
model.acceleration[:, node_index] = original_acceleration
# Inclined support is only applied in local X
- model.free_dofs[3 * node_index - 2] = false
+ model.free_dofs[3*node_index-2] = false
global_base = 3 * (node_index - 1) # Block index in global stiffness
model.global_transform[global_base+1:global_base+3, global_base+1:global_base+3] = bc.rotation_matrix
@@ -371,8 +371,8 @@ function find_point_in_mesh(
blk_id::Int,
tol::Float64
)
- node_indices, ξ, found = find_point_in_mesh(point,model.fom_model,blk_id,tol)
- return node_indices, ξ, found
+ node_indices, ξ, found = find_point_in_mesh(point, model.fom_model, blk_id, tol)
+ return node_indices, ξ, found
end
function apply_bc_detail(model::SolidMechanics, bc::SMContactSchwarzBC)
@@ -381,7 +381,7 @@ function apply_bc_detail(model::SolidMechanics, bc::SMContactSchwarzBC)
else
apply_sm_schwarz_contact_neumann(model, bc)
end
- if (bc.swap_bcs == true)
+ if (bc.swap_bcs == true)
bc.is_dirichlet = !bc.is_dirichlet
end
end
@@ -392,69 +392,69 @@ function apply_bc_detail(model::SolidMechanics, bc::CouplingSchwarzBoundaryCondi
else
apply_sm_schwarz_coupling_neumann(model, bc)
end
- if (bc.swap_bcs == true)
+ if (bc.swap_bcs == true)
bc.is_dirichlet = !bc.is_dirichlet
end
end
function apply_bc_detail(model::LinearOpInfRom, bc::CouplingSchwarzBoundaryCondition)
- if (typeof(bc.coupled_subsim.model) == SolidMechanics)
- ## Apply BC to the FOM vector
- apply_bc_detail(model.fom_model,bc)
-
- # populate our own BC vector
- bc_vector = zeros(3,length(bc.side_set_node_indices))
- for i ∈ 1:length(bc.side_set_node_indices)
- node_index = bc.side_set_node_indices[i]
- bc_vector[:,i] = model.fom_model.current[:, node_index] - model.fom_model.reference[:, node_index]
- end
- op_name = "B_"*bc.side_set_name
- bc_operator = model.opinf_rom[op_name]
- for i in 1:3
- model.reduced_boundary_forcing[:] += bc_operator[i,:,:] * bc_vector[i,:]
+ if (typeof(bc.coupled_subsim.model) == SolidMechanics)
+ ## Apply BC to the FOM vector
+ apply_bc_detail(model.fom_model, bc)
+
+ # populate our own BC vector
+ bc_vector = zeros(3, length(bc.side_set_node_indices))
+ for i ∈ 1:length(bc.side_set_node_indices)
+ node_index = bc.side_set_node_indices[i]
+ bc_vector[:, i] = model.fom_model.current[:, node_index] - model.fom_model.reference[:, node_index]
+ end
+ op_name = "B_" * bc.side_set_name
+ bc_operator = model.opinf_rom[op_name]
+ for i in 1:3
+ model.reduced_boundary_forcing[:] += bc_operator[i, :, :] * bc_vector[i, :]
+ end
+ else
+ throw("ROM-ROM coupling not supported yet")
end
- else
- throw("ROM-ROM coupling not supported yet")
- end
end
function apply_sm_schwarz_coupling_dirichlet(model::SolidMechanics, bc::CouplingSchwarzBoundaryCondition)
- if (typeof(bc.coupled_subsim.model) == SolidMechanics)
- for i ∈ 1:length(bc.side_set_node_indices)
- node_index = bc.side_set_node_indices[i]
- coupled_node_indices = bc.coupled_nodes_indices[i]
- N = bc.interpolation_function_values[i]
- elem_posn = bc.coupled_subsim.model.current[:, coupled_node_indices]
- elem_velo = bc.coupled_subsim.model.velocity[:, coupled_node_indices]
- elem_acce = bc.coupled_subsim.model.acceleration[:, coupled_node_indices]
- 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
- dof_index = [3 * node_index - 2, 3 * node_index - 1, 3 * node_index]
- model.free_dofs[dof_index] .= false
- end
- elseif (typeof(bc.coupled_subsim.model) == LinearOpInfRom)
- for i ∈ 1:length(bc.side_set_node_indices)
- node_index = bc.side_set_node_indices[i]
- coupled_node_indices = bc.coupled_nodes_indices[i]
- N = bc.interpolation_function_values[i]
- elem_posn = bc.coupled_subsim.model.fom_model.current[:, coupled_node_indices]
- elem_velo = bc.coupled_subsim.model.fom_model.velocity[:, coupled_node_indices]
- elem_acce = bc.coupled_subsim.model.fom_model.acceleration[:, coupled_node_indices]
- point_posn = elem_posn * N
- point_velo = elem_velo * N
- point_acce = elem_acce * N
- model.current[:, node_index] = point_posn
- model.velocity[:, node_index] = point_velo
- model.acceleration[:, node_index] = point_acce
- dof_index = [3 * node_index - 2, 3 * node_index - 1, 3 * node_index]
- model.free_dofs[dof_index] .= false
+ if (typeof(bc.coupled_subsim.model) == SolidMechanics)
+ for i ∈ 1:length(bc.side_set_node_indices)
+ node_index = bc.side_set_node_indices[i]
+ coupled_node_indices = bc.coupled_nodes_indices[i]
+ N = bc.interpolation_function_values[i]
+ elem_posn = bc.coupled_subsim.model.current[:, coupled_node_indices]
+ elem_velo = bc.coupled_subsim.model.velocity[:, coupled_node_indices]
+ elem_acce = bc.coupled_subsim.model.acceleration[:, coupled_node_indices]
+ 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
+ dof_index = [3 * node_index - 2, 3 * node_index - 1, 3 * node_index]
+ model.free_dofs[dof_index] .= false
+ end
+ elseif (typeof(bc.coupled_subsim.model) == LinearOpInfRom)
+ for i ∈ 1:length(bc.side_set_node_indices)
+ node_index = bc.side_set_node_indices[i]
+ coupled_node_indices = bc.coupled_nodes_indices[i]
+ N = bc.interpolation_function_values[i]
+ elem_posn = bc.coupled_subsim.model.fom_model.current[:, coupled_node_indices]
+ elem_velo = bc.coupled_subsim.model.fom_model.velocity[:, coupled_node_indices]
+ elem_acce = bc.coupled_subsim.model.fom_model.acceleration[:, coupled_node_indices]
+ point_posn = elem_posn * N
+ point_velo = elem_velo * N
+ point_acce = elem_acce * N
+ model.current[:, node_index] = point_posn
+ model.velocity[:, node_index] = point_velo
+ model.acceleration[:, node_index] = point_acce
+ dof_index = [3 * node_index - 2, 3 * node_index - 1, 3 * node_index]
+ model.free_dofs[dof_index] .= false
+ end
end
- end
end
function apply_sm_schwarz_coupling_neumann(model::SolidMechanics, bc::CouplingSchwarzBoundaryCondition)
@@ -549,9 +549,6 @@ function apply_bc(model::SolidMechanics, bc::SchwarzBoundaryCondition)
)
end
-
-
-
function apply_bc(model::LinearOpInfRom, bc::SchwarzBoundaryCondition)
global_sim = bc.coupled_subsim.params["global_simulation"]
schwarz_controller = global_sim.schwarz_controller
@@ -586,9 +583,9 @@ function apply_bc(model::LinearOpInfRom, bc::SchwarzBoundaryCondition)
interp_∂Ω_f =
same_step == true ? ∂Ω_f_hist[end] : interpolate(time_hist, ∂Ω_f_hist, time)
if typeof(bc.coupled_subsim.model) == SolidMechanics
- bc.coupled_subsim.model.internal_force = interp_∂Ω_f
+ bc.coupled_subsim.model.internal_force = interp_∂Ω_f
elseif typeof(bc.coupled_subsim.model) == LinearOpInfRom
- bc.coupled_subsim.model.fom_model.internal_force = interp_∂Ω_f
+ bc.coupled_subsim.model.fom_model.internal_force = interp_∂Ω_f
end
if typeof(bc) == SMContactSchwarzBC || typeof(bc) == SMNonOverlapSchwarzBC
@@ -668,7 +665,7 @@ function apply_sm_schwarz_contact_dirichlet(model::SolidMechanics, bc::SMContact
θ = asin(s)
m = w / s
rv = θ * m
- # Rotation is converted via the psuedo vector to rotation matrix
+ # Rotation is converted via the pseudo vector to rotation matrix
bc.rotation_matrix = MiniTensor.rt_from_rv(rv)
end
@@ -852,7 +849,7 @@ function create_bcs(params::Dict{String,Any})
inclined_support_nodes = Vector{Int64}()
for (bc_type, bc_type_params) ∈ bc_params
for bc_setting_params ∈ bc_type_params
- if bc_type == "Dirichlet"
+ if bc_type == "Dirichlet"
boundary_condition = SMDirichletBC(input_mesh, bc_setting_params)
push!(boundary_conditions, boundary_condition)
elseif bc_type == "Neumann"
@@ -914,7 +911,6 @@ function apply_bcs(model::LinearOpInfRom)
end
-
function apply_ics(params::Dict{String,Any}, model::SolidMechanics)
if haskey(params, "initial conditions") == false
return
@@ -954,31 +950,30 @@ end
function apply_ics(params::Dict{String,Any}, model::LinearOpInfRom)
- apply_ics(params,model.fom_model)
+ apply_ics(params, model.fom_model)
if haskey(params, "initial conditions") == false
return
end
- n_var,n_node,n_mode = model.basis.size
- n_var_fom,n_node_fom = size(model.fom_model.current)
+ n_var, n_node, n_mode = model.basis.size
+ n_var_fom, n_node_fom = size(model.fom_model.current)
# Make sure basis is the right size
- if n_var != n_var_fom || n_node != n_node_fom
- throw("Basis is wrong size")
+ if n_var != n_var_fom || n_node != n_node_fom
+ throw("Basis is wrong size")
end
# project onto basis
for k in 1:n_mode
- model.reduced_state[k] = 0.0
- for j in 1:n_node
- for n in 1:n_var
- model.reduced_state[k] += model.basis[n,j,k]*(model.fom_model.current[n,j] - model.fom_model.reference[n,j])
+ model.reduced_state[k] = 0.0
+ for j in 1:n_node
+ for n in 1:n_var
+ model.reduced_state[k] += model.basis[n, j, k] * (model.fom_model.current[n, j] - model.fom_model.reference[n, j])
+ end
end
- end
end
end
-
function pair_schwarz_bcs(sim::MultiDomainSimulation)
for subsim ∈ sim.subsims
model = subsim.model
diff --git a/src/model.jl b/src/model.jl
index 5b9f205..dffe26f 100644
--- a/src/model.jl
+++ b/src/model.jl
@@ -15,16 +15,8 @@ function LinearOpInfRom(params::Dict{String,Any})
opinf_model_file = params["model"]["model-file"]
opinf_model = NPZ.npzread(opinf_model_file)
basis = opinf_model["basis"]
- num_dofs_per_node,num_nodes_basis,reduced_dim = size(basis)
+ _, _, reduced_dim = size(basis)
num_dofs = reduced_dim
- #=
- coords = read_coordinates(fom_model.mesh)
- num_nodes = Exodus.num_nodes(fom_model.mesh.init)
- if (num_nodes_basis != num_nodes)
- throw("Basis size incompatible with mesh")
- end
- =#
-
time = 0.0
failed = false
null_vec = zeros(num_dofs)
@@ -45,12 +37,10 @@ function LinearOpInfRom(params::Dict{String,Any})
failed,
fom_model,
reference,
- false
+ false
)
end
-
-
function SolidMechanics(params::Dict{String,Any})
input_mesh = params["input_mesh"]
model_params = params["model"]
@@ -581,4 +571,4 @@ function evaluate(integrator::TimeIntegrator, model::SolidMechanics)
else
error("Unknown type of time integrator", typeof(integrator))
end
-end
+end
\ No newline at end of file
diff --git a/src/simulation.jl b/src/simulation.jl
index d450f12..166d219 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -42,7 +42,7 @@ function create_bcs(sim::SingleDomainSimulation)
end
end
sim.model.boundary_conditions = boundary_conditions
-
+
end
function create_bcs(sim::MultiDomainSimulation)
@@ -63,8 +63,8 @@ function SingleDomainSimulation(params::Dict{String,Any})
params["output_mesh"] = output_mesh
params["input_mesh"] = input_mesh
model = create_model(params)
- integrator = create_time_integrator(params,model)
- solver = create_solver(params,model)
+ integrator = create_time_integrator(params, model)
+ solver = create_solver(params, model)
failed = false
return SingleDomainSimulation(name, params, integrator, solver, model, failed)
end
diff --git a/src/solver.jl b/src/solver.jl
index 8b27392..d8fd5f4 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -17,7 +17,7 @@ function create_step(solver_params::Dict{String,Any})
end
end
-function HessianMinimizer(params::Dict{String,Any},model::Any)
+function HessianMinimizer(params::Dict{String,Any}, model::Any)
solver_params = params["solver"]
num_dof, = size(model.free_dofs)
minimum_iterations = solver_params["minimum iterations"]
@@ -167,11 +167,11 @@ function SteepestDescentStep(params::Dict{String,Any})
SteepestDescentStep(step_length)
end
-function create_solver(params::Dict{String,Any},model::Any)
+function create_solver(params::Dict{String,Any}, model::Any)
solver_params = params["solver"]
solver_name = solver_params["type"]
if solver_name == "Hessian minimizer"
- return HessianMinimizer(params,model)
+ return HessianMinimizer(params, model)
elseif solver_name == "explicit solver"
return ExplicitSolver(params)
elseif solver_name == "steepest descent"
@@ -233,27 +233,27 @@ function copy_solution_source_targets(
velocity = integrator.velocity
acceleration = integrator.acceleration
# Clean this up; maybe make a free dofs 2d array or move to a basis in matrix format
- for i = 1 : size(model.fom_model.current)[2]
- x_dof_index = 3 * (i - 1) + 1
- y_dof_index = 3 * (i - 1) + 2
- z_dof_index = 3 * (i - 1) + 3
- if model.fom_model.free_dofs[x_dof_index]
- model.fom_model.current[1,i] = model.basis[1,i,:]'displacement + model.fom_model.reference[1,i]
- model.fom_model.velocity[1,i] = model.basis[1,i,:]'velocity
- model.fom_model.acceleration[1,i] = model.basis[1,i,:]'acceleration
- end
+ for i = 1:size(model.fom_model.current)[2]
+ x_dof_index = 3 * (i - 1) + 1
+ y_dof_index = 3 * (i - 1) + 2
+ z_dof_index = 3 * (i - 1) + 3
+ if model.fom_model.free_dofs[x_dof_index]
+ model.fom_model.current[1, i] = model.basis[1, i, :]'displacement + model.fom_model.reference[1, i]
+ model.fom_model.velocity[1, i] = model.basis[1, i, :]'velocity
+ model.fom_model.acceleration[1, i] = model.basis[1, i, :]'acceleration
+ end
+
+ if model.fom_model.free_dofs[y_dof_index]
+ model.fom_model.current[2, i] = model.basis[2, i, :]'displacement + model.fom_model.reference[2, i]
+ model.fom_model.velocity[2, i] = model.basis[2, i, :]'velocity
+ model.fom_model.acceleration[2, i] = model.basis[2, i, :]'acceleration
+ end
- if model.fom_model.free_dofs[y_dof_index]
- model.fom_model.current[2,i] = model.basis[2,i,:]'displacement + model.fom_model.reference[2,i]
- model.fom_model.velocity[2,i] = model.basis[2,i,:]'velocity
- model.fom_model.acceleration[2,i] = model.basis[2,i,:]'acceleration
- end
-
- if model.fom_model.free_dofs[z_dof_index]
- model.fom_model.current[3,i] = model.basis[3,i,:]'displacement + model.fom_model.reference[3,i]
- model.fom_model.velocity[3,i] = model.basis[3,i,:]'velocity
- model.fom_model.acceleration[3,i] = model.basis[3,i,:]'acceleration
- end
+ if model.fom_model.free_dofs[z_dof_index]
+ model.fom_model.current[3, i] = model.basis[3, i, :]'displacement + model.fom_model.reference[3, i]
+ model.fom_model.velocity[3, i] = model.basis[3, i, :]'velocity
+ model.fom_model.acceleration[3, i] = model.basis[3, i, :]'acceleration
+ end
end
end
@@ -340,12 +340,12 @@ function copy_solution_source_targets(
nodal_displacement = model.current[:, node] - model.reference[:, node]
nodal_velocity = model.velocity[:, node]
nodal_acceleration = model.acceleration[:, node]
-
+
if model.inclined_support == true
- base = 3*(node-1) # Block index in global stiffness
+ base = 3 * (node - 1) # Block index in global stiffness
local_transform = model.global_transform[base+1:base+3, base+1:base+3]
nodal_displacement = local_transform * nodal_displacement
- nodal_velocity = local_transform * nodal_velocity
+ nodal_velocity = local_transform * nodal_velocity
nodal_acceleration = local_transform * nodal_acceleration
end
integrator.displacement[3*node-2:3*node] = nodal_displacement
@@ -370,11 +370,11 @@ function copy_solution_source_targets(
nodal_velocity = velocity[3*node-2:3*node]
nodal_acceleration = acceleration[3*node-2:3*node]
if model.inclined_support == true
- base = 3*(node-1) # Block index in global stiffness
+ base = 3 * (node - 1) # Block index in global stiffness
# Local (integrator) to global (model), use transpose
local_transform = model.global_transform[base+1:base+3, base+1:base+3]'
nodal_displacement = local_transform * nodal_displacement
- nodal_velocity = local_transform * nodal_velocity
+ nodal_velocity = local_transform * nodal_velocity
nodal_acceleration = local_transform * nodal_acceleration
end
model.current[:, node] = model.reference[:, node] + nodal_displacement
@@ -398,11 +398,11 @@ function copy_solution_source_targets(
nodal_velocity = velocity[3*node-2:3*node]
nodal_acceleration = acceleration[3*node-2:3*node]
if model.inclined_support == true
- base = 3*(node-1) # Block index in global stiffness
+ base = 3 * (node - 1) # Block index in global stiffness
# Local (integrator) to global (model), use transpose
local_transform = model.global_transform[base+1:base+3, base+1:base+3]'
nodal_displacement = local_transform * nodal_displacement
- nodal_velocity = local_transform * nodal_velocity
+ nodal_velocity = local_transform * nodal_velocity
nodal_acceleration = local_transform * nodal_acceleration
end
model.current[:, node] = model.reference[:, node] + nodal_displacement
@@ -422,11 +422,11 @@ function copy_solution_source_targets(
nodal_velocity = model.velocity[:, node]
nodal_acceleration = model.acceleration[:, node]
if model.inclined_support == true
- base = 3*(node-1) # Block index in global stiffness
+ base = 3 * (node - 1) # Block index in global stiffness
# Global (model) to local (integrator)
local_transform = model.global_transform[base+1:base+3, base+1:base+3]
nodal_displacement = local_transform * nodal_displacement
- nodal_velocity = local_transform * nodal_velocity
+ nodal_velocity = local_transform * nodal_velocity
nodal_acceleration = local_transform * nodal_acceleration
end
integrator.displacement[3*node-2:3*node] = nodal_displacement
@@ -439,7 +439,7 @@ end
### Move to model?
function evaluate(integrator::Newmark, solver::HessianMinimizer, model::LinearOpInfRom)
- beta = integrator.β
+ beta = integrator.β
gamma = integrator.γ
dt = integrator.time_step
@@ -452,13 +452,13 @@ function evaluate(integrator::Newmark, solver::HessianMinimizer, model::LinearOp
##M uddot + Ku = f
num_dof, = size(model.free_dofs)
- I = Matrix{Float64}(LinearAlgebra.I, num_dof,num_dof)
- LHS = I / (dt*dt*beta) + Matrix{Float64}(model.opinf_rom["K"])
- RHS = model.opinf_rom["f"] + model.reduced_boundary_forcing + 1.0/(dt*dt*beta).*integrator.disp_pre
+ I = Matrix{Float64}(LinearAlgebra.I, num_dof, num_dof)
+ LHS = I / (dt * dt * beta) + Matrix{Float64}(model.opinf_rom["K"])
+ RHS = model.opinf_rom["f"] + model.reduced_boundary_forcing + 1.0 / (dt * dt * beta) .* integrator.disp_pre
- residual = RHS - LHS * solver.solution
- solver.hessian[:,:] = LHS
- solver.gradient[:] = -residual
+ residual = RHS - LHS * solver.solution
+ solver.hessian[:, :] = LHS
+ solver.gradient[:] = -residual
end
@@ -476,7 +476,7 @@ function evaluate(integrator::QuasiStatic, solver::HessianMinimizer, model::Soli
else
solver.gradient = internal_force - external_force
solver.hessian = stiffness_matrix
- end
+ end
end
@@ -604,7 +604,7 @@ function compute_step(
solver::HessianMinimizer,
_::NewtonStep,
)
- return -solver.hessian\ solver.gradient
+ return -solver.hessian \ solver.gradient
end
diff --git a/src/time_integrator.jl b/src/time_integrator.jl
index 5d218a8..4f64bb4 100644
--- a/src/time_integrator.jl
+++ b/src/time_integrator.jl
@@ -67,7 +67,7 @@ function QuasiStatic(params::Dict{String,Any})
)
end
-function Newmark(params::Dict{String,Any},model::Any)
+function Newmark(params::Dict{String,Any}, model::Any)
integrator_params = params["time integrator"]
initial_time = integrator_params["initial time"]
final_time = integrator_params["final time"]
@@ -155,13 +155,13 @@ function CentralDifference(params::Dict{String,Any})
)
end
-function create_time_integrator(params::Dict{String,Any},model::Any)
+function create_time_integrator(params::Dict{String,Any}, model::Any)
integrator_params = params["time integrator"]
integrator_name = integrator_params["type"]
if integrator_name == "quasi static"
return QuasiStatic(params)
elseif integrator_name == "Newmark"
- return Newmark(params,model)
+ return Newmark(params, model)
elseif integrator_name == "central difference"
return CentralDifference(params)
else
@@ -199,8 +199,8 @@ function predict(integrator::Newmark, solver::Any, model::LinearOpInfRom)
dt = integrator.time_step
beta = integrator.β
gamma = integrator.γ
- integrator.disp_pre[:] = integrator.displacement[:] = integrator.displacement[:] + dt*integrator.velocity + 1.0/2.0*dt*dt*(1.0 - 2.0*beta)*integrator.acceleration
- integrator.velo_pre[:] = integrator.velocity[:] += dt*(1.0 - gamma)*integrator.acceleration
+ integrator.disp_pre[:] = integrator.displacement[:] = integrator.displacement[:] + dt * integrator.velocity + 1.0 / 2.0 * dt * dt * (1.0 - 2.0 * beta) * integrator.acceleration
+ integrator.velo_pre[:] = integrator.velocity[:] += dt * (1.0 - gamma) * integrator.acceleration
solver.solution[:] = integrator.displacement[:]
model.reduced_state[:] = integrator.displacement[:]
end
@@ -212,8 +212,8 @@ function correct(integrator::Newmark, solver::Any, model::LinearOpInfRom)
gamma = integrator.γ
integrator.displacement[:] = solver.solution[:]
- integrator.acceleration[:] = (solver.solution - integrator.disp_pre)/(beta*dt*dt)
- integrator.velocity[:] = integrator.velo_pre + dt*gamma*integrator.acceleration[:]
+ integrator.acceleration[:] = (solver.solution - integrator.disp_pre) / (beta * dt * dt)
+ integrator.velocity[:] = integrator.velo_pre + dt * gamma * integrator.acceleration[:]
model.reduced_state[:] = solver.solution[:]
@@ -350,7 +350,7 @@ function initialize_writing(
integrator::DynamicTimeIntegrator,
model::LinearOpInfRom,
)
- initialize_writing(params,integrator,model.fom_model)
+ initialize_writing(params, integrator, model.fom_model)
end
function initialize_writing(params::Dict{String,Any}, integrator::TimeIntegrator, _::SolidMechanics)
@@ -496,7 +496,7 @@ function write_step(params::Dict{String,Any}, integrator::Any, model::Any)
end
write_step_csv(integrator, model, sim_id)
if haskey(params, "CSV write sidesets") == true
- write_sideset_step_csv(params,integrator,model,sim_id)
+ write_sideset_step_csv(params, integrator, model, sim_id)
end
end
end
@@ -529,57 +529,57 @@ function write_step_csv(integrator::TimeIntegrator, model::SolidMechanics, sim_i
end
end
-function write_sideset_step_csv(params::Dict{String,Any},integrator::DynamicTimeIntegrator, model::SolidMechanics, sim_id::Integer)
- stop = integrator.stop
- index_string = "-" * string(stop, pad = 4)
- sim_id_string = string(sim_id, pad = 2) * "-"
- input_mesh = params["input_mesh"]
- bc_params = params["boundary conditions"]
-
- for bc ∈ model.boundary_conditions
- if (typeof(bc) == SMDirichletBC)
- node_set_name = bc.node_set_name
- offset = bc.offset
- if offset == 1
- offset_name = "x"
- end
- if offset == 2
- offset_name = "y"
- end
- if offset == 3
- offset_name = "z"
+function write_sideset_step_csv(params::Dict{String,Any}, integrator::DynamicTimeIntegrator, model::SolidMechanics, sim_id::Integer)
+ stop = integrator.stop
+ index_string = "-" * string(stop, pad=4)
+ sim_id_string = string(sim_id, pad=2) * "-"
+ input_mesh = params["input_mesh"]
+ bc_params = params["boundary conditions"]
+
+ for bc ∈ model.boundary_conditions
+ if (typeof(bc) == SMDirichletBC)
+ node_set_name = bc.node_set_name
+ offset = bc.offset
+ if offset == 1
+ offset_name = "x"
+ end
+ if offset == 2
+ offset_name = "y"
+ end
+ if offset == 3
+ offset_name = "z"
+ end
+ curr_filename = sim_id_string * node_set_name * "-" * offset_name * "-curr" * index_string * ".csv"
+ disp_filename = sim_id_string * node_set_name * "-" * offset_name * "-disp" * index_string * ".csv"
+ velo_filename = sim_id_string * node_set_name * "-" * offset_name * "-velo" * index_string * ".csv"
+ acce_filename = sim_id_string * node_set_name * "-" * offset_name * "-acce" * index_string * ".csv"
+ writedlm(curr_filename, model.current[bc.offset, bc.node_set_node_indices])
+ writedlm(velo_filename, model.velocity[bc.offset, bc.node_set_node_indices])
+ writedlm(acce_filename, model.acceleration[bc.offset, bc.node_set_node_indices])
+ writedlm(disp_filename, model.current[bc.offset, bc.node_set_node_indices] - model.reference[bc.offset, bc.node_set_node_indices])
+ elseif (typeof(bc) == SMOverlapSchwarzBC)
+ side_set_name = bc.side_set_name
+ curr_filename = sim_id_string * side_set_name * "-curr" * index_string * ".csv"
+ disp_filename = sim_id_string * side_set_name * "-disp" * index_string * ".csv"
+ velo_filename = sim_id_string * side_set_name * "-velo" * index_string * ".csv"
+ acce_filename = sim_id_string * side_set_name * "-acce" * index_string * ".csv"
+ writedlm_nodal_array(curr_filename, model.current[:, bc.side_set_node_indices])
+ writedlm_nodal_array(velo_filename, model.velocity[:, bc.side_set_node_indices])
+ writedlm_nodal_array(acce_filename, model.acceleration[:, bc.side_set_node_indices])
+ writedlm_nodal_array(disp_filename, model.current[:, bc.side_set_node_indices] - model.reference[:, bc.side_set_node_indices])
end
- curr_filename = sim_id_string * node_set_name * "-" * offset_name * "-curr" * index_string * ".csv"
- disp_filename = sim_id_string * node_set_name * "-" * offset_name * "-disp" * index_string * ".csv"
- velo_filename = sim_id_string * node_set_name * "-" * offset_name * "-velo" * index_string * ".csv"
- acce_filename = sim_id_string * node_set_name * "-" * offset_name * "-acce" * index_string * ".csv"
- writedlm(curr_filename, model.current[bc.offset,bc.node_set_node_indices])
- writedlm(velo_filename, model.velocity[bc.offset,bc.node_set_node_indices])
- writedlm(acce_filename, model.acceleration[bc.offset,bc.node_set_node_indices])
- writedlm(disp_filename, model.current[bc.offset,bc.node_set_node_indices] - model.reference[bc.offset,bc.node_set_node_indices])
- elseif (typeof(bc) == SMOverlapSchwarzBC)
- side_set_name = bc.side_set_name
- curr_filename = sim_id_string * side_set_name * "-curr" * index_string * ".csv"
- disp_filename = sim_id_string * side_set_name * "-disp" * index_string * ".csv"
- velo_filename = sim_id_string * side_set_name * "-velo" * index_string * ".csv"
- acce_filename = sim_id_string * side_set_name * "-acce" * index_string * ".csv"
- writedlm_nodal_array(curr_filename, model.current[:,bc.side_set_node_indices])
- writedlm_nodal_array(velo_filename, model.velocity[:,bc.side_set_node_indices])
- writedlm_nodal_array(acce_filename, model.acceleration[:,bc.side_set_node_indices])
- writedlm_nodal_array(disp_filename, model.current[:,bc.side_set_node_indices] - model.reference[:,bc.side_set_node_indices])
- end
- end
+ end
end
function write_step_csv(integrator::DynamicTimeIntegrator, model::OpInfModel, sim_id::Integer)
stop = integrator.stop
- index_string = "-" * string(stop, pad = 4)
- sim_id_string = string(sim_id, pad = 2) * "-"
+ index_string = "-" * string(stop, pad=4)
+ sim_id_string = string(sim_id, pad=2) * "-"
reduced_states_filename = sim_id_string * "reduced_states" * index_string * ".csv"
time_filename = sim_id_string * "time" * index_string * ".csv"
writedlm(reduced_states_filename, model.reduced_state)
- write_step_csv(integrator,model.fom_model,sim_id)
+ write_step_csv(integrator, model.fom_model, sim_id)
end
function write_step_exodus(
@@ -589,23 +589,23 @@ function write_step_exodus(
)
#Re-construct full state
reduced_state = model.reduced_state[:]
- for i = 1 : size(model.fom_model.current)[2]
- x_dof_index = 3 * (i - 1) + 1
- y_dof_index = 3 * (i - 1) + 2
- z_dof_index = 3 * (i - 1) + 3
- if model.fom_model.free_dofs[x_dof_index]
- model.fom_model.current[1,i] = model.basis[1,i,:]'reduced_state + model.fom_model.reference[1,i]
- end
-
- if model.fom_model.free_dofs[y_dof_index]
- model.fom_model.current[2,i] = model.basis[2,i,:]'reduced_state + model.fom_model.reference[2,i]
- end
-
- if model.fom_model.free_dofs[z_dof_index]
- model.fom_model.current[3,i] = model.basis[3,i,:]'reduced_state+ model.fom_model.reference[3,i]
- end
+ for i = 1:size(model.fom_model.current)[2]
+ x_dof_index = 3 * (i - 1) + 1
+ y_dof_index = 3 * (i - 1) + 2
+ z_dof_index = 3 * (i - 1) + 3
+ if model.fom_model.free_dofs[x_dof_index]
+ model.fom_model.current[1, i] = model.basis[1, i, :]'reduced_state + model.fom_model.reference[1, i]
+ end
+
+ if model.fom_model.free_dofs[y_dof_index]
+ model.fom_model.current[2, i] = model.basis[2, i, :]'reduced_state + model.fom_model.reference[2, i]
+ end
+
+ if model.fom_model.free_dofs[z_dof_index]
+ model.fom_model.current[3, i] = model.basis[3, i, :]'reduced_state + model.fom_model.reference[3, i]
+ end
end
- write_step_exodus(params,integrator,model.fom_model)
+ write_step_exodus(params, integrator, model.fom_model)
end