suanPan
MatrixModifier.hpp
Go to the documentation of this file.
1/*******************************************************************************
2 * Copyright (C) 2017-2024 Theodore Chang
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 ******************************************************************************/
30#ifndef MATRIXMODIFIER_H
31#define MATRIXMODIFIER_H
32
33#include <Element/Element.h>
34#include <suanPan.h>
35
36namespace suanpan {
37 namespace mass {
39 template<typename T> static void apply(Mat<T>&);
40 };
41
42 struct lumped_scale {
43 template<typename T> static void apply(Mat<T>&, unsigned);
44 };
45 } // namespace mass
46 namespace damping {
47 struct rayleigh {
48 template<typename T> static void apply(const shared_ptr<Element>&, T, T, T, T);
49 };
50
51 struct elemental {
52 template<typename T> static void apply(const shared_ptr<Element>&, T);
53 };
54 } // namespace damping
55} // namespace suanpan
56
57template<typename T> void suanpan::mass::lumped_simple::apply(Mat<T>& mass) { mass = diagmat(sum(mass)); }
58
59template<typename T> void suanpan::mass::lumped_scale::apply(Mat<T>& mass, const unsigned dim) {
60 Col<T> diag_mass(mass.n_rows, fill::zeros);
61
62 for(unsigned I = 0; I < dim; ++I) {
63 auto total_mass = 0.;
64 auto true_mass = 0.;
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));
68 }
69 if(fabs(true_mass) > 1E-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;
72 }
73 }
74
75 mass = diagmat(diag_mass);
76}
77
78template<typename T> void suanpan::damping::rayleigh::apply(const shared_ptr<Element>& element_obj, const T alpha, const T beta, const T zeta, const T eta) {
79 auto& ele_damping = access::rw(element_obj->get_trial_viscous());
80
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);
83
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();
88 }
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();
93 }
94
95 ele_damping = damping;
96
97 ele_force = damping * get_trial_velocity(element_obj.get());
98 }
99 else if(!ele_damping.is_empty()) ele_force = ele_damping * get_trial_velocity(element_obj.get());
100}
101
102template<typename T> void suanpan::damping::elemental::apply(const shared_ptr<Element>& element_obj, const T damping_ratio) {
103 const auto& t_stiffness = element_obj->get_current_stiffness();
104 const auto& t_mass = element_obj->get_current_mass();
105
106 if(t_stiffness.is_empty() || t_mass.is_empty()) return;
107
108 cx_vec eig_val;
109 cx_mat eig_vec;
110
111 auto num_mode = std::min(damping_ratio.n_elem, rank(t_mass));
112
113 if(!eig_pair(eig_val, eig_vec, t_stiffness, t_mass)) return;
114
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();
117
118 access::rw(element_obj->get_trial_viscous()) = damping;
119
120 access::rw(element_obj->get_trial_damping_force()) = damping * get_trial_velocity(element_obj.get());
121}
122
123#endif
124
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