libMesh
newmark_system.C
Go to the documentation of this file.
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
35 
36 
37 
39  const std::string & name_in,
40  const unsigned int number_in) :
41  LinearImplicitSystem (es, name_in, number_in),
42  _a_0 (1./(_default_alpha*_default_timestep*_default_timestep)),
43  _a_1 (_default_delta/(_default_alpha*_default_timestep)),
44  _a_2 (1./(_default_alpha*_default_timestep)),
45  _a_3 (1./(2.*_default_alpha)-1.),
46  _a_4 (_default_delta/_default_alpha -1.),
47  _a_5 (_default_timestep/2.*(_default_delta/_default_alpha-2.)),
48  _a_6 (_default_timestep*(1.-_default_delta)),
49  _a_7 (_default_delta*_default_timestep),
50  _finished_assemble (false)
51 
52 {
53  // default values of the newmark parameters
54  es.parameters.set<Real>("Newmark alpha") = _default_alpha;
55  es.parameters.set<Real>("Newmark delta") = _default_delta;
56 
57  // time step size.
58  // should be handled at a later stage through EquationSystems?
59  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  this->add_matrix ("stiffness");
69  this->add_matrix ("damping");
70  this->add_matrix ("mass");
71 
72  // load vector
73  this->add_vector ("force");
74 
75  // the displacement and the time derivatives
76  this->add_vector ("displacement");
77  this->add_vector ("velocity");
78  this->add_vector ("acceleration");
79 
80  // contributions to the rhs
81  this->add_vector ("rhs_m");
82  this->add_vector ("rhs_c");
83 
84  // results from the previous time step
85  this->add_vector ("old_solution");
86  this->add_vector ("old_acceleration");
87 }
88 
89 
90 
92 
93 
94 
96 {
97  // use parent clear this will also clear the
98  // matrices and vectors added in the constructor
100 
101  // Get a reference to the EquationSystems
102  EquationSystems & es =
103  this->get_equation_systems();
104 
105  // default values of the newmark parameters
106  es.parameters.set<Real>("Newmark alpha") = _default_alpha;
107  es.parameters.set<Real>("Newmark delta") = _default_delta;
108 
109  // time step size. should be handled at a later stage through EquationSystems?
110  es.parameters.set<Real>("Newmark time step") = _default_timestep;
111 
112  // set bool to false
113  _finished_assemble = false;
114 }
115 
116 
117 
119 {
120  libmesh_not_implemented();
121 }
122 
123 
124 
126 {
127  if (!_finished_assemble)
128  {
129  // prepare matrix with the help of the _dof_map,
130  // fill with sparsity pattern
132 
133  // compute the effective system matrix
134  this->compute_matrix();
135 
136  // apply initial conditions
137  this->initial_conditions();
138 
139  _finished_assemble = true;
140  }
141 }
142 
143 
145 {
146  // libmesh_assert(init_cond_fptr);
147 
148  // Log how long the user's matrix assembly code takes
149  LOG_SCOPE("initial_conditions ()", "NewmarkSystem");
150 
151  // Set all values to 0, then
152  // call the user-specified function for initial conditions.
153  this->get_vector("displacement").zero();
154  this->get_vector("velocity").zero();
155  this->get_vector("acceleration").zero();
156  this->user_initialization();
157 }
158 
159 
160 
162 {
163  // close the component matrices
164  this->get_matrix ("stiffness").close();
165  this->get_matrix ("mass" ).close();
166  this->get_matrix ("damping" ).close();
167 
168  // close & zero the system matrix
169  this->matrix->close (); this->matrix->zero();
170 
171  // add up the matrices
172  this->matrix->add (1., this->get_matrix ("stiffness"));
173  this->matrix->add (_a_0, this->get_matrix ("mass"));
174  this->matrix->add (_a_1, this->get_matrix ("damping"));
175 
176 }
177 
178 
179 
181 {
182  LOG_SCOPE("update_rhs ()", "NewmarkSystem");
183 
184  // zero the rhs-vector
185  NumericVector<Number> & the_rhs = *this->rhs;
186  the_rhs.zero();
187 
188  // get writable references to some vectors
189  NumericVector<Number> & rhs_m = this->get_vector("rhs_m");
190  NumericVector<Number> & rhs_c = this->get_vector("rhs_c");
191 
192  // zero the vectors for matrix-vector product
193  rhs_m.zero();
194  rhs_c.zero();
195 
196  // compute auxiliary vectors rhs_m and rhs_c
197  rhs_m.add(_a_0, this->get_vector("displacement"));
198  rhs_m.add(_a_2, this->get_vector("velocity"));
199  rhs_m.add(_a_3, this->get_vector("acceleration"));
200 
201  rhs_c.add(_a_1, this->get_vector("displacement"));
202  rhs_c.add(_a_4, this->get_vector("velocity"));
203  rhs_c.add(_a_5, this->get_vector("acceleration"));
204 
205  // compute rhs
206  NumericVector<Number> & force = this->get_vector("force");
207  force.close();
208  the_rhs.add(force);
209  the_rhs.add_vector(rhs_m, this->get_matrix("mass"));
210  the_rhs.add_vector(rhs_c, this->get_matrix("damping"));
211 }
212 
213 
214 
216 {
217  LOG_SCOPE("update_u_v_a ()", "NewmarkSystem");
218 
219  // get some references for convenience
220  const NumericVector<Number> & solu = *this->solution;
221 
222  NumericVector<Number> & disp_vec = this->get_vector("displacement");
223  NumericVector<Number> & vel_vec = this->get_vector("velocity");
224  NumericVector<Number> & acc_vec = this->get_vector("acceleration");
225  NumericVector<Number> & old_acc = this->get_vector("old_acceleration");
226  NumericVector<Number> & old_solu = this->get_vector("old_solution");
227 
228  // copy data
229  old_solu = disp_vec;
230  disp_vec = solu;
231  old_acc = acc_vec;
232 
233  // compute the new acceleration vector
234  acc_vec.scale(-_a_3);
235  acc_vec.add(_a_0, disp_vec);
236  acc_vec.add(-_a_0, old_solu);
237  acc_vec.add(-_a_2,vel_vec);
238 
239  // compute the new velocity vector
240  vel_vec.add(_a_6,old_acc);
241  vel_vec.add(_a_7,acc_vec);
242 }
243 
244 
245 
247  const Real alpha,
248  const Real delta)
249 {
250  libmesh_assert_not_equal_to (delta_T, 0.);
251 
252  // Get a reference to the EquationSystems
253  EquationSystems & es =
254  this->get_equation_systems();
255 
256  // the newmark parameters
257  es.parameters.set<Real>("Newmark alpha") = alpha;
258  es.parameters.set<Real>("Newmark delta") = delta;
259 
260  // time step size.
261  // should be handled at a later stage through EquationSystems?
262  es.parameters.set<Real>("Newmark time step") = delta_T;
263 
264  // the constants for time integration
265  _a_0 = 1./(alpha*delta_T*delta_T);
266  _a_1 = delta/(alpha*delta_T);
267  _a_2 = 1./(alpha*delta_T);
268  _a_3 = 1./(2.*alpha)-1.;
269  _a_4 = delta/alpha -1.;
270  _a_5 = delta_T/2.*(delta/alpha-2.);
271  _a_6 = delta_T*(1.-delta);
272  _a_7 = delta*delta_T;
273 }
274 
275 } // namespace libMesh
virtual void assemble() override
Prepares matrix and _dof_map for matrix assembly.
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...
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
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]...
const EquationSystems & get_equation_systems() const
Definition: system.h:721
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:2297
bool _finished_assemble
true if the matrix assembly is finished.
NumericVector< Number > & add_vector(std::string_view vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
Definition: system.C:768
The libMesh namespace provides an interface to certain functionality in the library.
virtual void assemble() override
Assemble the linear system.
virtual void zero()=0
Set all entries to zero.
virtual void add(const numeric_index_type i, const numeric_index_type j, const T value)=0
Add value to the element (i,j).
virtual void clear() override
Clear all the data structures associated with the system.
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.
virtual void zero()=0
Set all entries to 0.
virtual void reinit() override
Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be use...
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1593
static const Real _default_delta
Default Newmark delta.
void initial_conditions()
Apply initial conditions.
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
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.
Parameters parameters
Data structure holding arbitrary parameters.
virtual void clear() override
Clear all the data structures associated with the system.
virtual void add(const numeric_index_type i, const T value)=0
Adds value to the vector entry specified by i.
SparseMatrix< Number > & add_matrix(std::string_view mat_name, ParallelType type=PARALLEL, MatrixBuildType mat_build_type=MatrixBuildType::AUTOMATIC)
Adds the additional matrix mat_name to this system.
Definition: system.C:1010
const SparseMatrix< Number > & get_matrix(std::string_view mat_name) const
Definition: system.C:1124
const NumericVector< Number > & get_vector(std::string_view vec_name) const
Definition: system.C:943