libMesh
Functions
fem_system_ex3.C File Reference

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 58 of file fem_system_ex3.C.

References libMesh::DiffSolver::absolute_residual_tolerance, libMesh::BoundaryInfo::add_edge(), libMesh::BoundaryInfo::add_node(), libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), libMesh::MeshTools::Generation::build_cube(), libMesh::LibMeshInit::comm(), libMesh::NewmarkSolver::compute_initial_accel(), libMesh::default_solver_package(), libMesh::DifferentiableSystem::deltat, dim, libMesh::MeshBase::element_ptr_range(), libMesh::FIRST, libMesh::MeshBase::get_boundary_info(), libMesh::NewtonSolver::get_linear_solver(), libMesh::System::get_vector(), libMesh::BoundaryInfo::has_boundary_id(), libMesh::HEX8, libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::DiffSolver::initial_linear_tolerance, libMesh::INVALID_SOLVER_PACKAGE, libMesh::LAGRANGE, libMesh::DiffSolver::max_linear_iterations, libMesh::DiffSolver::max_nonlinear_iterations, mesh, libMesh::out, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::DiffSolver::quiet, libMesh::Real, libMesh::DiffSolver::relative_residual_tolerance, libMesh::DiffSolver::relative_step_tolerance, libMesh::LinearSolver< T >::set_solver_type(), side, libMesh::System::solution, libMesh::FEMSystem::solve(), libMesh::solver, libMesh::SPARSELU, libMesh::EulerSolver::theta, libMesh::Euler2Solver::theta, libMesh::System::time, libMesh::DifferentiableSystem::time_solver, libMesh::DiffSolver::verbose, and libMesh::ExodusII_IO::write_timestep().

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.
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
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
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.
This class implements a TimeSolver which does a single solve of the steady state problem.
Definition: steady_solver.h:47
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
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...