Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library 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 GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local includes
21 : #include "libmesh/newmark_system.h"
22 : #include "libmesh/equation_systems.h"
23 : #include "libmesh/sparse_matrix.h"
24 : #include "libmesh/libmesh_logging.h"
25 : #include "libmesh/numeric_vector.h"
26 :
27 : namespace libMesh
28 : {
29 :
30 : // ------------------------------------------------------------
31 : // NewmarkSystem static members
32 : const Real NewmarkSystem::_default_alpha = .25;
33 : const Real NewmarkSystem::_default_delta = .5;
34 : const Real NewmarkSystem::_default_timestep = 1.;
35 :
36 :
37 :
38 71 : NewmarkSystem::NewmarkSystem (EquationSystems & es,
39 : const std::string & name_in,
40 71 : const unsigned int number_in) :
41 : LinearImplicitSystem (es, name_in, number_in),
42 67 : _a_0 (1./(_default_alpha*_default_timestep*_default_timestep)),
43 67 : _a_1 (_default_delta/(_default_alpha*_default_timestep)),
44 67 : _a_2 (1./(_default_alpha*_default_timestep)),
45 67 : _a_3 (1./(2.*_default_alpha)-1.),
46 67 : _a_4 (_default_delta/_default_alpha -1.),
47 67 : _a_5 (_default_timestep/2.*(_default_delta/_default_alpha-2.)),
48 67 : _a_6 (_default_timestep*(1.-_default_delta)),
49 67 : _a_7 (_default_delta*_default_timestep),
50 71 : _finished_assemble (false)
51 :
52 : {
53 : // default values of the newmark parameters
54 71 : es.parameters.set<Real>("Newmark alpha") = _default_alpha;
55 71 : es.parameters.set<Real>("Newmark delta") = _default_delta;
56 :
57 : // time step size.
58 : // should be handled at a later stage through EquationSystems?
59 71 : es.parameters.set<Real>("Newmark time step") = _default_timestep;
60 :
61 : // add additional matrices and vectors that will be used in the
62 : // newmark algorithm to the data structure
63 : // functions LinearImplicitSystem::add_matrix and LinearImplicitSystem::add_vector
64 : // are used so we do not have to bother about initialization and
65 : // dof mapping
66 :
67 : // system matrices
68 71 : this->add_matrix ("stiffness");
69 71 : this->add_matrix ("damping");
70 71 : this->add_matrix ("mass");
71 :
72 : // load vector
73 71 : this->add_vector ("force");
74 :
75 : // the displacement and the time derivatives
76 71 : this->add_vector ("displacement");
77 71 : this->add_vector ("velocity");
78 71 : this->add_vector ("acceleration");
79 :
80 : // contributions to the rhs
81 71 : this->add_vector ("rhs_m");
82 71 : this->add_vector ("rhs_c");
83 :
84 : // results from the previous time step
85 71 : this->add_vector ("old_solution");
86 71 : this->add_vector ("old_acceleration");
87 71 : }
88 :
89 :
90 :
91 134 : NewmarkSystem::~NewmarkSystem () = default;
92 :
93 :
94 :
95 0 : void NewmarkSystem::clear ()
96 : {
97 : // use parent clear this will also clear the
98 : // matrices and vectors added in the constructor
99 0 : LinearImplicitSystem::clear();
100 :
101 : // Get a reference to the EquationSystems
102 : EquationSystems & es =
103 0 : this->get_equation_systems();
104 :
105 : // default values of the newmark parameters
106 0 : es.parameters.set<Real>("Newmark alpha") = _default_alpha;
107 0 : es.parameters.set<Real>("Newmark delta") = _default_delta;
108 :
109 : // time step size. should be handled at a later stage through EquationSystems?
110 0 : es.parameters.set<Real>("Newmark time step") = _default_timestep;
111 :
112 : // set bool to false
113 0 : _finished_assemble = false;
114 0 : }
115 :
116 :
117 :
118 0 : void NewmarkSystem::reinit ()
119 : {
120 0 : libmesh_not_implemented();
121 : }
122 :
123 :
124 :
125 21371 : void NewmarkSystem::assemble ()
126 : {
127 21371 : if (!_finished_assemble)
128 : {
129 : // prepare matrix with the help of the _dof_map,
130 : // fill with sparsity pattern
131 2 : LinearImplicitSystem::assemble();
132 :
133 : // compute the effective system matrix
134 71 : this->compute_matrix();
135 :
136 : // apply initial conditions
137 71 : this->initial_conditions();
138 :
139 71 : _finished_assemble = true;
140 : }
141 21371 : }
142 :
143 :
144 71 : void NewmarkSystem::initial_conditions ()
145 : {
146 : // libmesh_assert(init_cond_fptr);
147 :
148 : // Log how long the user's matrix assembly code takes
149 4 : LOG_SCOPE("initial_conditions ()", "NewmarkSystem");
150 :
151 : // Set all values to 0, then
152 : // call the user-specified function for initial conditions.
153 71 : this->get_vector("displacement").zero();
154 71 : this->get_vector("velocity").zero();
155 71 : this->get_vector("acceleration").zero();
156 71 : this->user_initialization();
157 71 : }
158 :
159 :
160 :
161 71 : void NewmarkSystem::compute_matrix ()
162 : {
163 : // close the component matrices
164 71 : this->get_matrix ("stiffness").close();
165 71 : this->get_matrix ("mass" ).close();
166 71 : this->get_matrix ("damping" ).close();
167 :
168 : // close & zero the system matrix
169 71 : this->matrix->close (); this->matrix->zero();
170 :
171 : // add up the matrices
172 71 : this->matrix->add (1., this->get_matrix ("stiffness"));
173 71 : this->matrix->add (_a_0, this->get_matrix ("mass"));
174 71 : this->matrix->add (_a_1, this->get_matrix ("damping"));
175 :
176 71 : }
177 :
178 :
179 :
180 21300 : void NewmarkSystem::update_rhs ()
181 : {
182 600 : LOG_SCOPE("update_rhs ()", "NewmarkSystem");
183 :
184 : // zero the rhs-vector
185 21300 : NumericVector<Number> & the_rhs = *this->rhs;
186 21300 : the_rhs.zero();
187 :
188 : // get writable references to some vectors
189 21300 : NumericVector<Number> & rhs_m = this->get_vector("rhs_m");
190 21300 : NumericVector<Number> & rhs_c = this->get_vector("rhs_c");
191 :
192 : // zero the vectors for matrix-vector product
193 21300 : rhs_m.zero();
194 21300 : rhs_c.zero();
195 :
196 : // compute auxiliary vectors rhs_m and rhs_c
197 21300 : rhs_m.add(_a_0, this->get_vector("displacement"));
198 21300 : rhs_m.add(_a_2, this->get_vector("velocity"));
199 21300 : rhs_m.add(_a_3, this->get_vector("acceleration"));
200 :
201 21300 : rhs_c.add(_a_1, this->get_vector("displacement"));
202 21300 : rhs_c.add(_a_4, this->get_vector("velocity"));
203 21300 : rhs_c.add(_a_5, this->get_vector("acceleration"));
204 :
205 : // compute rhs
206 21300 : NumericVector<Number> & force = this->get_vector("force");
207 21300 : force.close();
208 21300 : the_rhs.add(force);
209 21300 : the_rhs.add_vector(rhs_m, this->get_matrix("mass"));
210 21300 : the_rhs.add_vector(rhs_c, this->get_matrix("damping"));
211 21300 : }
212 :
213 :
214 :
215 21300 : void NewmarkSystem::update_u_v_a ()
216 : {
217 1200 : LOG_SCOPE("update_u_v_a ()", "NewmarkSystem");
218 :
219 : // get some references for convenience
220 600 : const NumericVector<Number> & solu = *this->solution;
221 :
222 21300 : NumericVector<Number> & disp_vec = this->get_vector("displacement");
223 21300 : NumericVector<Number> & vel_vec = this->get_vector("velocity");
224 21300 : NumericVector<Number> & acc_vec = this->get_vector("acceleration");
225 21300 : NumericVector<Number> & old_acc = this->get_vector("old_acceleration");
226 21300 : NumericVector<Number> & old_solu = this->get_vector("old_solution");
227 :
228 : // copy data
229 21300 : old_solu = disp_vec;
230 21300 : disp_vec = solu;
231 21300 : old_acc = acc_vec;
232 :
233 : // compute the new acceleration vector
234 21300 : acc_vec.scale(-_a_3);
235 21300 : acc_vec.add(_a_0, disp_vec);
236 21300 : acc_vec.add(-_a_0, old_solu);
237 21300 : acc_vec.add(-_a_2,vel_vec);
238 :
239 : // compute the new velocity vector
240 21300 : vel_vec.add(_a_6,old_acc);
241 21300 : vel_vec.add(_a_7,acc_vec);
242 21300 : }
243 :
244 :
245 :
246 71 : void NewmarkSystem::set_newmark_parameters (const Real delta_T,
247 : const Real alpha,
248 : const Real delta)
249 : {
250 2 : libmesh_assert_not_equal_to (delta_T, 0.);
251 :
252 : // Get a reference to the EquationSystems
253 : EquationSystems & es =
254 4 : this->get_equation_systems();
255 :
256 : // the newmark parameters
257 71 : es.parameters.set<Real>("Newmark alpha") = alpha;
258 71 : es.parameters.set<Real>("Newmark delta") = delta;
259 :
260 : // time step size.
261 : // should be handled at a later stage through EquationSystems?
262 71 : es.parameters.set<Real>("Newmark time step") = delta_T;
263 :
264 : // the constants for time integration
265 71 : _a_0 = 1./(alpha*delta_T*delta_T);
266 71 : _a_1 = delta/(alpha*delta_T);
267 71 : _a_2 = 1./(alpha*delta_T);
268 71 : _a_3 = 1./(2.*alpha)-1.;
269 71 : _a_4 = delta/alpha -1.;
270 71 : _a_5 = delta_T/2.*(delta/alpha-2.);
271 71 : _a_6 = delta_T*(1.-delta);
272 71 : _a_7 = delta*delta_T;
273 71 : }
274 :
275 : } // namespace libMesh
|
