Skip to content

Commit

Permalink
Use new string element in the bow simulation
Browse files Browse the repository at this point in the history
  • Loading branch information
@stfnp
stfnp committed Oct 16, 2024
1 parent 95ba4fa commit 0b934fe
Show file tree
Hide file tree
Showing 6 changed files with 87 additions and 95 deletions.
8 changes: 4 additions & 4 deletions solver/bows/users/2r77c5r2.bow
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,13 +42,13 @@
0.1,
0.01
],
[
1.0,
0.006
],
[
0.2,
0.005
],
[
1.0,
0.006
]
],
"name": "ydin",
Expand Down
92 changes: 38 additions & 54 deletions solver/src/bow/simulation.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ use iter_num_tools::lin_space;
use itertools::Itertools;
use nalgebra::vector;
use soa_rs::Soa;
use crate::fem::solvers::eigen::{Mode, natural_frequencies, natural_frequencies_from_matrices};
use crate::fem::solvers::eigen::{Mode, natural_frequencies};
use crate::fem::solvers::statics::StaticSolver;
use crate::fem::system::element::Element;
use crate::fem::system::node::Node;
Expand Down Expand Up @@ -39,7 +39,7 @@ pub struct Simulation<'a> {
limb_nodes: Vec<Node>,
limb_elements: Vec<usize>,

string_node: Node,
string_nodes: Vec<Node>,
string_element: usize,

// Arrow mass element positioned at the string center
Expand Down Expand Up @@ -96,8 +96,6 @@ impl<'a> Simulation<'a> {
for &e in &limb_elements {
system.element_mut::<BeamElement>(e).set_damping(alpha);
}

println!("{}, {}", -1.0, modes[0].omega);
}

// Layer setup data
Expand All @@ -115,46 +113,44 @@ impl<'a> Simulation<'a> {
let limb_width = s_eval.iter().map(|&s| { section.width(s) }).collect();
let limb_height = s_eval.iter().map(|&s| { section.height(s) }).collect();

// Limb tip and string center positions
let x_tip = u_nodes.last().unwrap()[0];
let y_tip = u_nodes.last().unwrap()[1];
let y_str = -input.dimensions.brace_height;
// String center node that is fixed in the case of no string.
// The rest of the string nodes come from the limb.
let string_center = system.create_node(&vector![0.0, -input.dimensions.brace_height, 0.0], &[false, string, false]);
let mut string_nodes = vec![string_center];
string_nodes.extend_from_slice(&limb_nodes);

// Add string center node and make it fixed in the case of no string
let string_node = if string {
system.create_node(&vector![0.0, y_str, 0.0], &[false, true, false])
}
else {
system.create_node(&vector![0.0, y_str, 0.0], &[false, false, false])
};
// Arrow mass element is placed at the string center
let arrow_element = system.add_element(&[string_nodes[0]], MassElement::new(0.5*input.masses.arrow));

// The string element only gets non-zero parameters if the string option is true
let l = f64::hypot(x_tip, y_str - y_tip);
let EA = if string { (input.string.n_strands as f64)*input.string.strand_stiffness } else { 0.0 };
let ρA = if string { (input.string.n_strands as f64)*input.string.strand_density } else { 0.0 };
let string_element = StringElement::bar(EA, 0.0, l); // Damping is determined later when the length of the string is known
let string_element = system.add_element(&[*limb_nodes.last().unwrap(), string_node], string_element);

// Arrow mass and other remaining additional masses
let arrow_element = system.add_element(&[string_node], MassElement::new(0.5*input.masses.arrow));
system.add_element(&[string_node], MassElement::new(0.5*input.masses.string_center + 1.0/3.0*ρA*l));
system.add_element(&[*limb_nodes.last().unwrap()], MassElement::new(input.masses.string_tip + 2.0/3.0*ρA*l));
let offsets = vec![0.0; limb_nodes.len() + 1];
let string_element = StringElement::new(EA, 0.0, 1.0, offsets); // Damping is determined later when the length of the string is known
let string_element = system.add_element(&string_nodes, string_element);

// Evaluate the string element so that the actual string length is computed
// Then set the initial length to the current length so that the string is tension-free.
system.eval_element(string_element);
let element = system.element_mut::<StringElement>(string_element);
let l0 = element.get_current_length();
element.set_initial_length(l0);

// Finish the simulation info object
let simulation = Self {
input,
limb_nodes,
limb_elements,
string_node,
string_nodes,
string_element,
arrow_element,
arrow_separation: None,
};

// If string is to be initialized, perform bracing simulation
if string {
let l0 = system.element_ref::<StringElement>(string_element).get_initial_length();
system.add_force(string_node.y(), move |_t| { -1.0 });
system.add_force(simulation.string_nodes[0].y(), move |_t| { -1.0 });

// Initial values for the string factor, the slope and the step size for iterating on the string factor
let mut factor1 = 1.0;
Expand All @@ -174,7 +170,7 @@ impl<'a> Simulation<'a> {
system.element_mut::<StringElement>(simulation.string_element).set_initial_length(factor*l0);

let mut solver = StaticSolver::new(&mut system, statics::Settings::default()); // TODO: Don't construct new solver in each iteration
let result = solver.equilibrium_displacement_controlled(simulation.string_node.y(), -input.dimensions.brace_height);
let result = solver.equilibrium_displacement_controlled(simulation.string_nodes[0].y(), -input.dimensions.brace_height);
let slope = simulation.get_string_slope(&system);

(slope, result)
Expand Down Expand Up @@ -215,9 +211,13 @@ impl<'a> Simulation<'a> {

// After the string length is known, we can calculate the viscosity that is required
// for achieving the prescribed string damping ratio
let l = system.element_ref::<StringElement>(simulation.string_element).get_initial_length();
let ηA = 4.0*l/PI*f64::sqrt(ρA*EA)*input.damping.damping_ratio_string;
let l0 = system.element_ref::<StringElement>(simulation.string_element).get_initial_length();
let ηA = 4.0*l0/PI*f64::sqrt(ρA*EA)*input.damping.damping_ratio_string;
system.element_mut::<StringElement>(simulation.string_element).set_linear_damping(ηA);

// Base + additional masses of the string
system.add_element(&[simulation.string_nodes[0]], MassElement::new(0.5*input.masses.string_center + 1.0/3.0*ρA*l0));
system.add_element(&[*simulation.limb_nodes.last().unwrap()], MassElement::new(input.masses.string_tip + 2.0/3.0*ρA*l0));
}

let common = Common {
Expand Down Expand Up @@ -249,7 +249,7 @@ impl<'a> Simulation<'a> {
let mut states = Soa::<State>::new();
let mut solver = StaticSolver::new(&mut system, statics::Settings::default());

solver.equilibrium_path_displacement_controlled(simulation.string_node.y(), -model.dimensions.draw_length, model.settings.n_draw_steps, &mut |system, eval| {
solver.equilibrium_path_displacement_controlled(simulation.string_nodes[0].y(), -model.dimensions.draw_length, model.settings.n_draw_steps, &mut |system, eval| {
let state = simulation.get_bow_state(&system, SystemEval::Static(&eval));
let progress = 100.0*(state.draw_length - model.dimensions.brace_height)/(model.dimensions.draw_length - model.dimensions.brace_height);
states.push(state);
Expand Down Expand Up @@ -314,10 +314,6 @@ impl<'a> Simulation<'a> {
// Push bow state into the final results
states.push(state);

let modes = natural_frequencies_from_matrices(eval.get_mass_matrix(), eval.get_damping_matrix(), eval.get_stiffness_matrix()).unwrap();
println!("{}, {}", system.get_time(), modes[0].omega);


// Continue the simulation as long as the arrow is not separated
// and the timeout has not yet been reached
return simulation.arrow_separation.is_none() && system.get_time() < t_max;
Expand All @@ -338,9 +334,6 @@ impl<'a> Simulation<'a> {
let state = simulation.get_bow_state(&system, SystemEval::Dynamic(&eval));
states.push(state);

let modes = natural_frequencies_from_matrices(eval.get_mass_matrix(), eval.get_damping_matrix(), eval.get_stiffness_matrix()).unwrap();
println!("{}, {}", system.get_time(), modes[0].omega);

// Continue as long as the end time is not reached
return system.get_time() < t_end;
}).map_err(|e| ModelError::SimulationDynamicSolutionFailed(e))?;
Expand Down Expand Up @@ -416,10 +409,10 @@ impl<'a> Simulation<'a> {
// TODO: Let mutation, write functions for intermediate results
fn get_bow_state(&self, system: &System, eval: SystemEval) -> State {
let time = system.get_time();
let draw_length = -system.get_displacement(self.string_node.y());
let draw_length = -system.get_displacement(self.string_nodes[0].y());
let draw_force = match eval {
SystemEval::Static(eval) => -2.0*eval.get_scaled_external_force(self.string_node.y()),
SystemEval::Dynamic(eval) => -2.0*eval.get_external_force(self.string_node.y())
SystemEval::Static(eval) => -2.0*eval.get_scaled_external_force(self.string_nodes[0].y()),
SystemEval::Dynamic(eval) => -2.0*eval.get_external_force(self.string_nodes[0].y())
};

let string_force = system.element_ref::<StringElement>(self.string_element).normal_force();
Expand All @@ -439,26 +432,17 @@ impl<'a> Simulation<'a> {
(
match eval {
SystemEval::Static(_) => 0.0,
SystemEval::Dynamic(eval) => eval.get_acceleration(self.string_node.y())
SystemEval::Dynamic(eval) => eval.get_acceleration(self.string_nodes[0].y())
},
system.get_velocity(self.string_node.y()),
system.get_displacement(self.string_node.y()),
system.get_velocity(self.string_nodes[0].y()),
system.get_displacement(self.string_nodes[0].y()),
)
};

// String kinematics

let limb_tip = self.limb_nodes.last().unwrap();

let string_pos = vec![
[ system.get_displacement(limb_tip.x()), system.get_displacement(limb_tip.y()) ],
[ system.get_displacement(self.string_node.x()), system.get_displacement(self.string_node.y()) ],
];

let string_vel = vec![
[ system.get_velocity(limb_tip.x()), system.get_velocity(limb_tip.y()) ],
[ system.get_velocity(self.string_node.x()), system.get_velocity(self.string_node.y()) ],
];
let string_pos = system.element_ref::<StringElement>(self.string_element).contact_points().map(|p| p.into()).collect_vec();
let string_vel = vec![[0.0, 0.0]; string_pos.len()]; // TODO

/*
let acc_string = self.string_node.map(|node| vec![
Expand Down Expand Up @@ -534,7 +518,7 @@ impl<'a> Simulation<'a> {
fn get_string_slope(&self, system: &System) -> f64 {
let x_tip = system.get_displacement(self.limb_nodes.last().unwrap().x());
let y_tip = system.get_displacement(self.limb_nodes.last().unwrap().y());
let y_str = system.get_displacement(self.string_node.y());
let y_str = system.get_displacement(self.string_nodes[0].y());
(y_tip - y_str)/x_tip
}
}
Expand Down
3 changes: 0 additions & 3 deletions solver/src/fem/elements/beam/linear.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -192,9 +192,6 @@ mod tests {
let C_sec = section.C(0.0);
let K_ref = LinearBeamSegmentFEM::element_stiffness_matrix(C_sec[(0, 0)], C_sec[(1, 1)], C_sec[(2, 2)], curve.length(), alpha);

println!("Error K_num = {}", (segment_num.K - K_ref).component_div(&K_ref));
println!("Error K_fem = {}", (segment_num.K - K_ref).component_div(&K_ref));

// Check stiffness matrices
assert_relative_eq!(segment_num.K, K_ref, max_relative=1e-9);
assert_relative_eq!(segment_fem.K, K_ref, max_relative=1e-6);
Expand Down
7 changes: 7 additions & 0 deletions solver/src/fem/elements/string.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,10 @@ impl StringElement {
Self::bar(k*l0, d*l0, l0)
}

pub fn get_current_length(&self) -> f64 {
self.lt
}

pub fn get_initial_length(&self) -> f64 {
self.l0
}
Expand Down Expand Up @@ -83,6 +87,9 @@ impl StringElement {
self.Nv
}

pub fn contact_points(&self) -> impl Iterator<Item=SVector<f64, 2>> + '_ {
self.indices.iter().map(|&i| self.points[i])
}
}

impl Element for StringElement {
Expand Down
10 changes: 10 additions & 0 deletions solver/src/fem/system/system.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -296,6 +296,16 @@ impl System {
eval
}

// Evaluates the element with the current system state
// Only necessary in special occasions, usually the solvers do this anyway.
pub fn eval_element(&mut self, index: usize) {
let (dofs, element) = &mut self.elements[index];
let u_view = PositionView::new(&self.u, dofs);
let v_view = VelocityView::new(&self.v, dofs);

element.set_state_and_evaluate(&u_view, &v_view, None, None, None);
}

// TODO: Unify with other eval functions
pub fn eval_statics(&mut self, eval: &mut StaticEval) {
// Set to zero in case of previous values
Expand Down
Loading

0 comments on commit 0b934fe

Please sign in to comment.