libMesh
newmark_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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"
24 #include "libmesh/equation_systems.h"
25 #include "libmesh/sparse_matrix.h"
26 #include "libmesh/libmesh_logging.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_not_implemented();
131 }
132 
133 
134 
136 {
137  if (!_finished_assemble)
138  {
139  // prepare matrix with the help of the _dof_map,
140  // fill with sparsity pattern
142 
143  // compute the effective system matrix
144  this->compute_matrix();
145 
146  // apply initial conditions
147  this->initial_conditions();
148 
149  _finished_assemble = true;
150  }
151 }
152 
153 
155 {
156  // libmesh_assert(init_cond_fptr);
157 
158  // Log how long the user's matrix assembly code takes
159  LOG_SCOPE("initial_conditions ()", "NewmarkSystem");
160 
161  // Set all values to 0, then
162  // call the user-specified function for initial conditions.
163  this->get_vector("displacement").zero();
164  this->get_vector("velocity").zero();
165  this->get_vector("acceleration").zero();
166  this->user_initialization();
167 }
168 
169 
170 
172 {
173  // close the component matrices
174  this->get_matrix ("stiffness").close();
175  this->get_matrix ("mass" ).close();
176  this->get_matrix ("damping" ).close();
177 
178  // close & zero the system matrix
179  this->matrix->close (); this->matrix->zero();
180 
181  // add up the matrices
182  this->matrix->add (1., this->get_matrix ("stiffness"));
183  this->matrix->add (_a_0, this->get_matrix ("mass"));
184  this->matrix->add (_a_1, this->get_matrix ("damping"));
185 
186 }
187 
188 
189 
191 {
192  LOG_SCOPE("update_rhs ()", "NewmarkSystem");
193 
194  // zero the rhs-vector
195  NumericVector<Number> & the_rhs = *this->rhs;
196  the_rhs.zero();
197 
198  // get writable references to some vectors
199  NumericVector<Number> & rhs_m = this->get_vector("rhs_m");
200  NumericVector<Number> & rhs_c = this->get_vector("rhs_c");
201 
202  // zero the vectors for matrix-vector product
203  rhs_m.zero();
204  rhs_c.zero();
205 
206  // compute auxiliary vectors rhs_m and rhs_c
207  rhs_m.add(_a_0, this->get_vector("displacement"));
208  rhs_m.add(_a_2, this->get_vector("velocity"));
209  rhs_m.add(_a_3, this->get_vector("acceleration"));
210 
211  rhs_c.add(_a_1, this->get_vector("displacement"));
212  rhs_c.add(_a_4, this->get_vector("velocity"));
213  rhs_c.add(_a_5, this->get_vector("acceleration"));
214 
215  // compute rhs
216  the_rhs.add(this->get_vector("force"));
217  the_rhs.add_vector(rhs_m, this->get_matrix("mass"));
218  the_rhs.add_vector(rhs_c, this->get_matrix("damping"));
219 }
220 
221 
222 
224 {
225  LOG_SCOPE("update_u_v_a ()", "NewmarkSystem");
226 
227  // get some references for convenience
228  const NumericVector<Number> & solu = *this->solution;
229 
230  NumericVector<Number> & disp_vec = this->get_vector("displacement");
231  NumericVector<Number> & vel_vec = this->get_vector("velocity");
232  NumericVector<Number> & acc_vec = this->get_vector("acceleration");
233  NumericVector<Number> & old_acc = this->get_vector("old_acceleration");
234  NumericVector<Number> & old_solu = this->get_vector("old_solution");
235 
236  // copy data
237  old_solu = disp_vec;
238  disp_vec = solu;
239  old_acc = acc_vec;
240 
241  // compute the new acceleration vector
242  acc_vec.scale(-_a_3);
243  acc_vec.add(_a_0, disp_vec);
244  acc_vec.add(-_a_0, old_solu);
245  acc_vec.add(-_a_2,vel_vec);
246 
247  // compute the new velocity vector
248  vel_vec.add(_a_6,old_acc);
249  vel_vec.add(_a_7,acc_vec);
250 }
251 
252 
253 
255  const Real alpha,
256  const Real delta)
257 {
258  libmesh_assert_not_equal_to (delta_T, 0.);
259 
260  // Get a reference to the EquationSystems
261  EquationSystems & es =
262  this->get_equation_systems();
263 
264  // the newmark parameters
265  es.parameters.set<Real>("Newmark alpha") = alpha;
266  es.parameters.set<Real>("Newmark delta") = delta;
267 
268  // time step size.
269  // should be handled at a later stage through EquationSystems?
270  es.parameters.set<Real>("Newmark time step") = delta_T;
271 
272  // the constants for time integration
273  _a_0 = 1./(alpha*delta_T*delta_T);
274  _a_1 = delta/(alpha*delta_T);
275  _a_2 = 1./(alpha*delta_T);
276  _a_3 = 1./(2.*alpha)-1.;
277  _a_4 = delta/alpha -1.;
278  _a_5 = delta_T/2.*(delta/alpha-2.);
279  _a_6 = delta_T*(1.-delta);
280  _a_7 = delta*delta_T;
281 }
282 
283 } // namespace libMesh
This is the EquationSystems class.
void set_newmark_parameters(const Real delta_T=_default_timestep, const Real alpha=_default_alpha, const Real delta=_default_delta)
Set the time step size and the newmark parameter alpha and delta and calculate the constant parameter...
~NewmarkSystem()
Destructor.
virtual void clear() libmesh_override
Clear all the data structures associated with the system.
This class provides a specific system class.
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
SparseMatrix< Number > & add_matrix(const std::string &mat_name)
Adds the additional matrix mat_name to this system.
NumericVector< Number > * rhs
The system matrix.
virtual void user_initialization()
Calls user&#39;s attached initialization function, or is overridden by the user in derived classes...
Definition: system.C:1936
bool _finished_assemble
true if the matrix assembly is finished.
The libMesh namespace provides an interface to certain functionality in the library.
virtual void zero()=0
Set all entries to zero.
virtual void clear() libmesh_override
Clear all the data structures associated with the system.
virtual void add(const numeric_index_type i, const numeric_index_type j, const T value)=0
Add value to the element (i,j).
void update_rhs()
Update the rhs.
Real _a_0
Constants used for the time integration.
virtual void scale(const T factor)=0
Scale each element of the vector by the given factor.
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
Definition: system.C:681
virtual void zero()=0
Set all entries to 0.
virtual void assemble() libmesh_override
Prepares matrix and _dof_map for matrix assembly.
static const Real _default_delta
Default Newmark delta.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:794
void initial_conditions()
Apply initial conditions.
virtual void reinit() libmesh_override
Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be use...
const SparseMatrix< Number > & get_matrix(const std::string &mat_name) const
void update_u_v_a()
Update displacement, velocity and acceleration.
virtual void close()=0
Calls the SparseMatrix&#39;s internal assembly routines, ensuring that the values are consistent across p...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
T & set(const std::string &)
Definition: parameters.h:469
SparseMatrix< Number > * matrix
The system matrix.
static const Real _default_timestep
Default Newmark time step.
NewmarkSystem(EquationSystems &es, const std::string &name, const unsigned int number)
Constructor.
void compute_matrix()
Compute the global matrix by adding up scaled mass damping and stiffness matrix.
static const Real _default_alpha
Default Newmark alpha.
const EquationSystems & get_equation_systems() const
Definition: system.h:712
Parameters parameters
Data structure holding arbitrary parameters.
virtual void assemble() libmesh_override
Assemble the linear system.
virtual void add(const numeric_index_type i, const T value)=0
Adds value to each entry of the vector.