newmark_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 // C++ includes
21 
22 // Local includes
23 #include "libmesh/newmark_system.h"
25 #include "libmesh/sparse_matrix.h"
27 #include "libmesh/numeric_vector.h"
28 
29 namespace libMesh
30 {
31 
32 
33 
34 
35 // ------------------------------------------------------------
36 // NewmarkSystem static members
40 
41 
42 
43 // ------------------------------------------------------------
44 // NewmarkSystem implementation
46  const std::string& name_in,
47  const unsigned int number_in) :
48  LinearImplicitSystem (es, name_in, number_in),
49  _a_0 (1./(_default_alpha*_default_timestep*_default_timestep)),
50  _a_1 (_default_delta/(_default_alpha*_default_timestep)),
51  _a_2 (1./(_default_alpha*_default_timestep)),
52  _a_3 (1./(2.*_default_alpha)-1.),
53  _a_4 (_default_delta/_default_alpha -1.),
54  _a_5 (_default_timestep/2.*(_default_delta/_default_alpha-2.)),
55  _a_6 (_default_timestep*(1.-_default_delta)),
56  _a_7 (_default_delta*_default_timestep),
57  _finished_assemble (false)
58 
59 {
60  // default values of the newmark parameters
61  es.parameters.set<Real>("Newmark alpha") = _default_alpha;
62  es.parameters.set<Real>("Newmark delta") = _default_delta;
63 
64  // time step size.
65  // should be handled at a later stage through EquationSystems?
66  es.parameters.set<Real>("Newmark time step") = _default_timestep;
67 
68  // add additional matrices and vectors that will be used in the
69  // newmark algorithm to the data structure
70  // functions LinearImplicitSystem::add_matrix and LinearImplicitSystem::add_vector
71  // are used so we do not have to bother about initialization and
72  // dof mapping
73 
74  // system matrices
75  this->add_matrix ("stiffness");
76  this->add_matrix ("damping");
77  this->add_matrix ("mass");
78 
79  // load vector
80  this->add_vector ("force");
81 
82  // the displacement and the time derivatives
83  this->add_vector ("displacement");
84  this->add_vector ("velocity");
85  this->add_vector ("acceleration");
86 
87  // contributions to the rhs
88  this->add_vector ("rhs_m");
89  this->add_vector ("rhs_c");
90 
91  // results from the previous time step
92  this->add_vector ("old_solution");
93  this->add_vector ("old_acceleration");
94 }
95 
96 
97 
99 {
100  this->clear();
101 }
102 
103 
104 
106 {
107  // use parent clear this will also clear the
108  // matrices and vectors added in the constructor
110 
111  // Get a reference to the EquationSystems
112  EquationSystems& es =
113  this->get_equation_systems();
114 
115  // default values of the newmark parameters
116  es.parameters.set<Real>("Newmark alpha") = _default_alpha;
117  es.parameters.set<Real>("Newmark delta") = _default_delta;
118 
119  // time step size. should be handled at a later stage through EquationSystems?
120  es.parameters.set<Real>("Newmark time step") = _default_timestep;
121 
122  // set bool to false
123  _finished_assemble = false;
124 }
125 
126 
127 
129 {
130  libmesh_error();
131 
132  // initialize parent data
134 }
135 
136 
137 
139 {
140  if (!_finished_assemble)
141  {
142  // prepare matrix with the help of the _dof_map,
143  // fill with sparsity pattern
145 
146  // compute the effective system matrix
147  this->compute_matrix();
148 
149  // apply initial conditions
150  this->initial_conditions();
151 
152  _finished_assemble = true;
153  }
154 }
155 
156 
158 {
159  // libmesh_assert(init_cond_fptr);
160 
161  // Log how long the user's matrix assembly code takes
162  START_LOG("initial_conditions ()", "NewmarkSystem");
163 
164  // Set all values to 0, then
165  // call the user-specified function for initial conditions.
166  this->get_vector("displacement").zero();
167  this->get_vector("velocity").zero();
168  this->get_vector("acceleration").zero();
169  this->user_initialization();
170 
171  // Stop logging the user code
172  STOP_LOG("initial_conditions ()", "NewmarkSystem");
173 }
174 
175 
176 
178 {
179  // close the component matrices
180  this->get_matrix ("stiffness").close();
181  this->get_matrix ("mass" ).close();
182  this->get_matrix ("damping" ).close();
183 
184  // close & zero the system matrix
185  this->matrix->close (); this->matrix->zero();
186 
187  // add up the matrices
188  this->matrix->add (1., this->get_matrix ("stiffness"));
189  this->matrix->add (_a_0, this->get_matrix ("mass"));
190  this->matrix->add (_a_1, this->get_matrix ("damping"));
191 
192 }
193 
194 
195 
197 {
198  START_LOG("update_rhs ()", "NewmarkSystem");
199 
200  // zero the rhs-vector
201  NumericVector<Number>& the_rhs = *this->rhs;
202  the_rhs.zero();
203 
204  // get writable references to some vectors
205  NumericVector<Number>& rhs_m = this->get_vector("rhs_m");
206  NumericVector<Number>& rhs_c = this->get_vector("rhs_c");
207 
208 
209  // zero the vectors for matrix-vector product
210  rhs_m.zero();
211  rhs_c.zero();
212 
213  // compute auxiliary vectors rhs_m and rhs_c
214  rhs_m.add(_a_0, this->get_vector("displacement"));
215  rhs_m.add(_a_2, this->get_vector("velocity"));
216  rhs_m.add(_a_3, this->get_vector("acceleration"));
217 
218  rhs_c.add(_a_1, this->get_vector("displacement"));
219  rhs_c.add(_a_4, this->get_vector("velocity"));
220  rhs_c.add(_a_5, this->get_vector("acceleration"));
221 
222  // compute rhs
223  the_rhs.add(this->get_vector("force"));
224  the_rhs.add_vector(rhs_m, this->get_matrix("mass"));
225  the_rhs.add_vector(rhs_c, this->get_matrix("damping"));
226 
227  STOP_LOG("update_rhs ()", "NewmarkSystem");
228 }
229 
230 
231 
233 {
234  START_LOG("update_u_v_a ()", "NewmarkSystem");
235 
236  // get some references for convenience
237  const NumericVector<Number>& solu = *this->solution;
238 
239  NumericVector<Number>& disp_vec = this->get_vector("displacement");
240  NumericVector<Number>& vel_vec = this->get_vector("velocity");
241  NumericVector<Number>& acc_vec = this->get_vector("acceleration");
242  NumericVector<Number>& old_acc = this->get_vector("old_acceleration");
243  NumericVector<Number>& old_solu = this->get_vector("old_solution");
244 
245  // copy data
246  old_solu = disp_vec;
247  disp_vec = solu;
248  old_acc = acc_vec;
249 
250  // compute the new acceleration vector
251  acc_vec.scale(-_a_3);
252  acc_vec.add(_a_0, disp_vec);
253  acc_vec.add(-_a_0, old_solu);
254  acc_vec.add(-_a_2,vel_vec);
255 
256  // compute the new velocity vector
257  vel_vec.add(_a_6,old_acc);
258  vel_vec.add(_a_7,acc_vec);
259 
260  STOP_LOG("update_u_v_a ()", "NewmarkSystem");
261 }
262 
263 
264 
266  const Real alpha,
267  const Real delta)
268 {
269  libmesh_assert_not_equal_to (delta_T, 0.);
270 
271  // Get a reference to the EquationSystems
272  EquationSystems& es =
273  this->get_equation_systems();
274 
275  // the newmark parameters
276  es.parameters.set<Real>("Newmark alpha") = alpha;
277  es.parameters.set<Real>("Newmark delta") = delta;
278 
279  // time step size.
280  // should be handled at a later stage through EquationSystems?
281  es.parameters.set<Real>("Newmark time step") = delta_T;
282 
283  // the constants for time integration
284  _a_0 = 1./(alpha*delta_T*delta_T);
285  _a_1 = delta/(alpha*delta_T);
286  _a_2 = 1./(alpha*delta_T);
287  _a_3 = 1./(2.*alpha)-1.;
288  _a_4 = delta/alpha -1.;
289  _a_5 = delta_T/2.*(delta/alpha-2.);
290  _a_6 = delta_T*(1.-delta);
291  _a_7 = delta*delta_T;
292 }
293 
294 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo