30#ifndef MATRIXMODIFIER_H
31#define MATRIXMODIFIER_H
39 template<
typename T>
static void apply(Mat<T>&);
43 template<
typename T>
static void apply(Mat<T>&,
unsigned);
48 template<
typename T>
static void apply(
const shared_ptr<Element>&,
T,
T,
T,
T);
52 template<
typename T>
static void apply(
const shared_ptr<Element>&,
T);
60 Col<T> diag_mass(mass.n_rows, fill::zeros);
62 for(
unsigned I = 0; I < dim; ++I) {
65 for(
auto J = I; J < diag_mass.n_elem; J += dim) {
66 true_mass += mass(J, J);
67 total_mass += sum(mass.row(J));
69 if(fabs(true_mass) > 1
E-14) {
70 auto factor = total_mass / true_mass;
71 for(
auto J = I; J < diag_mass.n_elem; J += dim) diag_mass(J) = mass(J, J) * factor;
75 mass = diagmat(diag_mass);
79 auto& ele_damping = access::rw(element_obj->get_trial_viscous());
81 if(
auto& ele_force = access::rw(element_obj->get_trial_damping_force()); element_obj->if_update_viscous()) {
82 mat damping(element_obj->get_total_number(), element_obj->get_total_number(), fill::zeros);
84 if(0. != alpha && !element_obj->get_current_mass().is_empty()) damping += alpha * element_obj->get_current_mass();
85 if(0. != beta && !element_obj->get_current_stiffness().is_empty()) {
86 damping += beta * element_obj->get_current_stiffness();
87 if(element_obj->is_nlgeom() && !element_obj->get_current_geometry().is_empty()) damping += beta * element_obj->get_current_geometry();
89 if(0. != zeta && !element_obj->get_initial_stiffness().is_empty()) damping += zeta * element_obj->get_initial_stiffness();
90 if(0. != eta && !element_obj->get_trial_stiffness().is_empty()) {
91 damping += eta * element_obj->get_trial_stiffness();
92 if(element_obj->is_nlgeom() && !element_obj->get_trial_geometry().is_empty()) damping += eta * element_obj->get_trial_geometry();
95 ele_damping = damping;
99 else if(!ele_damping.is_empty()) ele_force = ele_damping *
get_trial_velocity(element_obj.get());
103 const auto& t_stiffness = element_obj->get_current_stiffness();
104 const auto& t_mass = element_obj->get_current_mass();
106 if(t_stiffness.is_empty() || t_mass.is_empty())
return;
111 auto num_mode = std::min(damping_ratio.n_elem, rank(t_mass));
113 if(!eig_pair(eig_val, eig_vec, t_stiffness, t_mass))
return;
115 const mat theta = t_mass * abs(eig_vec.head_cols(num_mode));
116 const mat damping = theta * diagmat(2. * sqrt(abs(eig_val.head(num_mode)) * damping_ratio)) * theta.t();
118 access::rw(element_obj->get_trial_viscous()) = damping;
120 access::rw(element_obj->get_trial_damping_force()) = damping *
get_trial_velocity(element_obj.get());
vec get_trial_velocity(const shared_ptr< Node > &)
Definition: NodeHelper.hpp:38
static void apply(Mat< T > &, unsigned)
Definition: MatrixModifier.hpp:59
static void apply(Mat< T > &)
Definition: MatrixModifier.hpp:57
static void apply(const shared_ptr< Element > &, T)
Definition: MatrixModifier.hpp:102
static void apply(const shared_ptr< Element > &, T, T, T, T)
Definition: MatrixModifier.hpp:78
Definition: MatrixModifier.hpp:36
Definition: MatrixModifier.hpp:51
Definition: MatrixModifier.hpp:47
Definition: MatrixModifier.hpp:42
Definition: MatrixModifier.hpp:38