libMesh
fem_system_ex3.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 // <h1>FEMSystem Example 3 - Unsteady Linear Elasticity with
21 // FEMSystem</h1>
22 // \author Paul Bauman
23 // \date 2015
24 //
25 // This example shows how to solve the three-dimensional transient
26 // linear elasticity equations using the DifferentiableSystem class framework.
27 // This is just Systems of Equations Example 6 recast.
28 
29 // C++ includes
30 #include <iomanip>
31 
32 // Basic include files
33 #include "libmesh/equation_systems.h"
34 #include "libmesh/error_vector.h"
35 #include "libmesh/getpot.h"
36 #include "libmesh/exodusII_io.h"
37 #include "libmesh/kelly_error_estimator.h"
38 #include "libmesh/mesh.h"
39 #include "libmesh/mesh_generation.h"
40 
41 // The systems and solvers we may use
42 #include "elasticity_system.h"
43 #include "libmesh/diff_solver.h"
44 #include "libmesh/newmark_solver.h"
45 #include "libmesh/steady_solver.h"
46 #include "libmesh/euler_solver.h"
47 #include "libmesh/euler2_solver.h"
48 #include "libmesh/elem.h"
49 #include "libmesh/newton_solver.h"
50 #include "libmesh/eigen_sparse_linear_solver.h"
51 
52 #define x_scaling 1.3
53 
54 // Bring in everything from the libMesh namespace
55 using namespace libMesh;
56 
57 // The main program.
58 int main (int argc, char ** argv)
59 {
60  // Initialize libMesh.
61  LibMeshInit init (argc, argv);
62 
63  // This example requires a linear solver package.
64  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
65  "--enable-petsc, --enable-trilinos, or --enable-eigen");
66 
67  // Parse the input file
68  GetPot infile("fem_system_ex3.in");
69 
70  // Override input file arguments from the command line
71  infile.parse_command_line(argc, argv);
72 
73  // Read in parameters from the input file
74  const Real deltat = infile("deltat", 0.25);
75  unsigned int n_timesteps = infile("n_timesteps", 1);
76 
77 #ifdef LIBMESH_HAVE_EXODUS_API
78  const unsigned int write_interval = infile("write_interval", 1);
79 #endif
80 
81  // Initialize the cantilever mesh
82  const unsigned int dim = 3;
83 
84  // Make sure libMesh was compiled for 3D
85  libmesh_example_requires(dim == LIBMESH_DIM, "3D support");
86 
87  // Create a 3D mesh distributed across the default MPI communicator.
88  Mesh mesh(init.comm(), dim);
90  32,
91  8,
92  4,
93  0., 1.*x_scaling,
94  0., 0.3,
95  0., 0.1,
96  HEX8);
97 
98 
99  // Print information about the mesh to the screen.
100  mesh.print_info();
101 
102  // Let's add some node and edge boundary conditions.
103  // Each processor should know about each boundary condition it can
104  // see, so we loop over all elements, not just local elements.
105  for (const auto & elem : mesh.element_ptr_range())
106  {
107  unsigned int
108  side_max_x = 0, side_min_y = 0,
109  side_max_y = 0, side_max_z = 0;
110  bool
111  found_side_max_x = false, found_side_max_y = false,
112  found_side_min_y = false, found_side_max_z = false;
113  for (auto side : elem->side_index_range())
114  {
115  if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_X))
116  {
117  side_max_x = side;
118  found_side_max_x = true;
119  }
120 
121  if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MIN_Y))
122  {
123  side_min_y = side;
124  found_side_min_y = true;
125  }
126 
127  if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_Y))
128  {
129  side_max_y = side;
130  found_side_max_y = true;
131  }
132 
133  if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_Z))
134  {
135  side_max_z = side;
136  found_side_max_z = true;
137  }
138  }
139 
140  // If elem has sides on boundaries
141  // BOUNDARY_ID_MAX_X, BOUNDARY_ID_MAX_Y, BOUNDARY_ID_MAX_Z
142  // then let's set a node boundary condition
143  if (found_side_max_x && found_side_max_y && found_side_max_z)
144  for (auto n : elem->node_index_range())
145  if (elem->is_node_on_side(n, side_max_x) &&
146  elem->is_node_on_side(n, side_max_y) &&
147  elem->is_node_on_side(n, side_max_z))
148  mesh.get_boundary_info().add_node(elem->node_ptr(n), NODE_BOUNDARY_ID);
149 
150  // If elem has sides on boundaries
151  // BOUNDARY_ID_MAX_X and BOUNDARY_ID_MIN_Y
152  // then let's set an edge boundary condition
153  if (found_side_max_x && found_side_min_y)
154  for (auto e : elem->edge_index_range())
155  if (elem->is_edge_on_side(e, side_max_x) &&
156  elem->is_edge_on_side(e, side_min_y))
157  mesh.get_boundary_info().add_edge(elem, e, EDGE_BOUNDARY_ID);
158  }
159 
160  // Create an equation systems object.
161  EquationSystems equation_systems (mesh);
162 
163  // Declare the system "Navier-Stokes" and its variables.
164  ElasticitySystem & system =
165  equation_systems.add_system<ElasticitySystem> ("Linear Elasticity");
166 
167  // Solve this as a time-dependent or steady system
168  std::string time_solver = infile("time_solver","DIE!");
169 
170  ExplicitSystem * v_system;
171  ExplicitSystem * a_system;
172 
173  if( time_solver == std::string("newmark") )
174  {
175  // Create ExplicitSystem to help output velocity
176  v_system = &equation_systems.add_system<ExplicitSystem> ("Velocity");
177  v_system->add_variable("u_vel", FIRST, LAGRANGE);
178  v_system->add_variable("v_vel", FIRST, LAGRANGE);
179  v_system->add_variable("w_vel", FIRST, LAGRANGE);
180 
181  // Create ExplicitSystem to help output acceleration
182  a_system = &equation_systems.add_system<ExplicitSystem> ("Acceleration");
183  a_system->add_variable("u_accel", FIRST, LAGRANGE);
184  a_system->add_variable("v_accel", FIRST, LAGRANGE);
185  a_system->add_variable("w_accel", FIRST, LAGRANGE);
186  }
187 
188  if( time_solver == std::string("newmark"))
189  system.time_solver.reset(new NewmarkSolver(system));
190 
191  else if( time_solver == std::string("euler") )
192  {
193  system.time_solver.reset(new EulerSolver(system));
194  EulerSolver & euler_solver = cast_ref<EulerSolver &>(*(system.time_solver.get()));
195  euler_solver.theta = infile("theta", 1.0);
196  }
197 
198  else if( time_solver == std::string("euler2") )
199  {
200  system.time_solver.reset(new Euler2Solver(system));
201  Euler2Solver & euler_solver = cast_ref<Euler2Solver &>(*(system.time_solver.get()));
202  euler_solver.theta = infile("theta", 1.0);
203  }
204 
205  else if( time_solver == std::string("steady"))
206  {
207  system.time_solver.reset(new SteadySolver(system));
208  libmesh_assert_equal_to (n_timesteps, 1);
209  }
210  else
211  libmesh_error_msg(std::string("ERROR: invalid time_solver ")+time_solver);
212 
213  // Initialize the system
214  equation_systems.init ();
215 
216  // Set the time stepping options
217  system.deltat = deltat;
218 
219  // And the nonlinear solver options
220  DiffSolver & solver = *(system.time_solver->diff_solver().get());
221  solver.quiet = infile("solver_quiet", true);
222  solver.verbose = !solver.quiet;
223  solver.max_nonlinear_iterations = infile("max_nonlinear_iterations", 15);
224  solver.relative_step_tolerance = infile("relative_step_tolerance", 1.e-3);
225  solver.relative_residual_tolerance = infile("relative_residual_tolerance", 0.0);
226  solver.absolute_residual_tolerance = infile("absolute_residual_tolerance", 0.0);
227 
228  // And the linear solver options
229  solver.max_linear_iterations = infile("max_linear_iterations", 50000);
230  solver.initial_linear_tolerance = infile("initial_linear_tolerance", 1.e-3);
231 
232  // Print information about the system to the screen.
233  equation_systems.print_info();
234 
235  // If we're using EulerSolver or Euler2Solver, and we're using EigenSparseLinearSolver,
236  // then we need to reset the EigenSparseLinearSolver to use SPARSELU because BICGSTAB
237  // chokes on the Jacobian matrix we give it and Eigen's GMRES currently doesn't work.
238  NewtonSolver * newton_solver = dynamic_cast<NewtonSolver *>( &solver );
239  if( newton_solver &&
240  (time_solver == std::string("euler") || time_solver == std::string("euler2") ) )
241  {
242  LinearSolver<Number> & linear_solver = newton_solver->get_linear_solver();
243 #ifdef LIBMESH_HAVE_EIGEN_SPARSE
244  EigenSparseLinearSolver<Number> * eigen_linear_solver =
245  dynamic_cast<EigenSparseLinearSolver<Number> *>(&linear_solver);
246 
247  if( eigen_linear_solver )
248  eigen_linear_solver->set_solver_type(SPARSELU);
249 #endif
250  }
251 
252  if( time_solver == std::string("newmark") )
253  {
254  NewmarkSolver * newmark = cast_ptr<NewmarkSolver*>(system.time_solver.get());
255  newmark->compute_initial_accel();
256 
257  // Copy over initial velocity and acceleration for output.
258  // Note we can do this because of the matching variables/FE spaces
259  *(v_system->solution) = system.get_vector("_old_solution_rate");
260  *(a_system->solution) = system.get_vector("_old_solution_accel");
261  }
262 
263 #ifdef LIBMESH_HAVE_EXODUS_API
264  // Output initial state
265  {
266  std::ostringstream file_name;
267 
268  // We write the file in the ExodusII format.
269  file_name << std::string("out.")+time_solver+std::string(".e-s.")
270  << std::setw(3)
271  << std::setfill('0')
272  << std::right
273  << 0;
274 
275  ExodusII_IO(mesh).write_timestep(file_name.str(),
276  equation_systems,
277  1, // This number indicates how many time steps
278  // are being written to the file
279  system.time);
280  }
281 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
282 
283  // Now we begin the timestep loop to compute the time-accurate
284  // solution of the equations.
285  for (unsigned int t_step=0; t_step != n_timesteps; ++t_step)
286  {
287  // A pretty update message
288  libMesh::out << "\n\nSolving time step "
289  << t_step
290  << ", time = "
291  << system.time
292  << std::endl;
293 
294  system.solve();
295 
296  // Advance to the next timestep in a transient problem
297  system.time_solver->advance_timestep();
298 
299  // Copy over updated velocity and acceleration for output.
300  // Note we can do this because of the matching variables/FE spaces
301  if( time_solver == std::string("newmark") )
302  {
303  *(v_system->solution) = system.get_vector("_old_solution_rate");
304  *(a_system->solution) = system.get_vector("_old_solution_accel");
305  }
306 
307 #ifdef LIBMESH_HAVE_EXODUS_API
308  // Write out this timestep if we're requested to
309  if ((t_step+1)%write_interval == 0)
310  {
311  std::ostringstream file_name;
312 
313  // We write the file in the ExodusII format.
314  file_name << std::string("out.")+time_solver+std::string(".e-s.")
315  << std::setw(3)
316  << std::setfill('0')
317  << std::right
318  << t_step+1;
319 
320  ExodusII_IO(mesh).write_timestep(file_name.str(),
321  equation_systems,
322  1, // This number indicates how many time steps
323  // are being written to the file
324  system.time);
325  }
326 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
327  }
328 
329  // All done.
330  return 0;
331 }
Real time
For time-dependent problems, this is the time t at the beginning of the current timestep.
Definition: system.h:1545
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
UniquePtr< TimeSolver > time_solver
A pointer to the solver object we&#39;re going to use.
Definition: diff_system.h:221
LinearSolver< Number > & get_linear_solver()
Definition: newton_solver.h:81
This is the EquationSystems class.
int main(int argc, char **argv)
void write_timestep(const std::string &fname, const EquationSystems &es, const int timestep, const Real time)
Writes out the solution at a specific timestep.
Definition: exodusII_io.C:773
bool quiet
The DiffSolver should not print anything to libMesh::out unless quiet is set to false; default is tru...
Definition: diff_solver.h:162
unsigned int dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
Real theta
The value for the theta method to employ: 1.0 corresponds to backwards Euler, 0.0 corresponds to forw...
Definition: euler_solver.h:100
unsigned int max_nonlinear_iterations
The DiffSolver should exit in failure if max_nonlinear_iterations is exceeded and continue_after_max_...
Definition: diff_solver.h:156
Real absolute_residual_tolerance
The DiffSolver should exit after the residual is reduced to either less than absolute_residual_tolera...
Definition: diff_solver.h:191
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
unsigned short int side
Definition: xdr_io.C:49
MeshBase & mesh
This class provides an interface to Eigen iterative solvers that is compatible with the libMesh Linea...
This class defines a Newmark time integrator for second order (in time) DifferentiableSystems.
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
The libMesh namespace provides an interface to certain functionality in the library.
Real initial_linear_tolerance
Any required linear solves will at first be done with this tolerance; the DiffSolver may tighten the ...
Definition: diff_solver.h:210
This class defines a solver which uses the default libMesh linear solver in a quasiNewton method to h...
Definition: newton_solver.h:46
virtual void compute_initial_accel()
This method uses the specified initial displacement and velocity to compute the initial acceleration ...
PetscDiffSolver & solver
SolverPackage default_solver_package()
Definition: libmesh.C:995
This is a generic class that defines a solver to handle ImplicitSystem classes, including NonlinearIm...
Definition: diff_solver.h:70
virtual void solve() libmesh_override
Invokes the solver associated with the system.
Definition: fem_system.C:1056
void add_node(const Node *node, const boundary_id_type id)
Add Node node with boundary id id to the boundary information data structures.
This class defines a theta-method (defaulting to Backward Euler with theta = 1.0) solver to handle ti...
Definition: euler2_solver.h:48
virtual SimpleRange< element_iterator > element_ptr_range()=0
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
Real deltat
For time-dependent problems, this is the amount delta t to advance the solution in time...
Definition: diff_system.h:248
unsigned int max_linear_iterations
Each linear solver step should exit after max_linear_iterations is exceeded.
Definition: diff_solver.h:148
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
This class implements a TimeSolver which does a single solve of the steady state problem.
Definition: steady_solver.h:47
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
Real theta
The value for the theta method to employ: 1.0 corresponds to backwards Euler, 0.0 corresponds to forw...
void set_solver_type(const SolverType st)
Sets the type of solver to use.
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:794
This class defines a theta-method Euler (defaulting to Backward Euler with theta = 1...
Definition: euler_solver.h:43
bool verbose
The DiffSolver may print a lot more to libMesh::out if verbose is set to true; default is false...
Definition: diff_solver.h:168
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
OStreamProxy out
bool has_boundary_id(const Node *const node, const boundary_id_type id) const
virtual void init()
Initialize all the systems.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
void build_cube(UnstructuredMesh &mesh, const unsigned int nx=0, const unsigned int ny=0, const unsigned int nz=0, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const Real zmin=0., const Real zmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
Builds a (elements) cube.
The ExplicitSystem provides only "right hand side" storage, which should be sufficient for solving mo...
Real relative_residual_tolerance
Definition: diff_solver.h:192
void add_edge(const dof_id_type elem, const unsigned short int edge, const boundary_id_type id)
Add edge edge of element number elem with boundary id id to the boundary information data structure...