libMesh
fem_system_ex2.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 // <h1>FEMSystem Example 2</h1>
20 // \author Robert Weidlich
21 // \date 2012
22 
23 #include "libmesh/equation_systems.h"
24 #include "libmesh/getpot.h"
25 #include "libmesh/libmesh.h"
26 #include "libmesh/libmesh_logging.h"
27 #include "libmesh/mesh.h"
28 #include "libmesh/mesh_generation.h"
29 #include "libmesh/numeric_vector.h"
30 #include "libmesh/string_to_enum.h"
31 #include "libmesh/time_solver.h"
32 #include "libmesh/transient_system.h"
33 #include "libmesh/vtk_io.h"
34 
35 #include <cstdio>
36 #include <ctime>
37 #include <fstream>
38 #include <iomanip>
39 #include <iostream>
40 #include <sstream>
41 #include <string>
42 #include <unistd.h>
43 
44 using namespace libMesh;
45 
46 #include "solid_system.h"
47 
48 void setup(EquationSystems & systems,
49  Mesh & mesh,
50  GetPot & args)
51 {
52  const unsigned int dim = mesh.mesh_dimension();
53  // We currently invert tensors with the assumption that they're 3x3
54  libmesh_assert (dim == 3);
55 
56  // Generating Mesh
57  ElemType eltype = Utility::string_to_enum<ElemType>(args("mesh/generation/element_type", "hex8"));
58  int nx = args("mesh/generation/num_elem", 4, 0);
59  int ny = args("mesh/generation/num_elem", 4, 1);
60  int nz = dim > 2 ? args("mesh/generation/num_elem", 4, 2) : 0;
61  double origx = args("mesh/generation/origin", -1.0, 0);
62  double origy = args("mesh/generation/origin", -1.0, 1);
63  double origz = args("mesh/generation/origin", 0.0, 2);
64  double sizex = args("mesh/generation/size", 2.0, 0);
65  double sizey = args("mesh/generation/size", 2.0, 1);
66  double sizez = args("mesh/generation/size", 2.0, 2);
67  MeshTools::Generation::build_cube(mesh, nx, ny, nz,
68  origx, origx+sizex, origy, origy+sizey, origz, origz+sizez, eltype);
69 
70  // Creating Systems
71  SolidSystem & imms = systems.add_system<SolidSystem> ("solid");
72  imms.args = args;
73 
74  // Build up auxiliary system
75  ExplicitSystem & aux_sys = systems.add_system<TransientExplicitSystem>("auxiliary");
76 
77  // Initialize the system
78  systems.parameters.set<unsigned int>("phase") = 0;
79  systems.init();
80 
81  imms.save_initial_mesh();
82 
83  // Fill global solution vector from local ones
84  aux_sys.reinit();
85 }
86 
87 
88 
89 void run_timestepping(EquationSystems & systems, GetPot & args)
90 {
91  TransientExplicitSystem & aux_system = systems.get_system<TransientExplicitSystem>("auxiliary");
92 
93  SolidSystem & solid_system = systems.get_system<SolidSystem>("solid");
94 
95  UniquePtr<VTKIO> io = UniquePtr<VTKIO>(new VTKIO(systems.get_mesh()));
96 
97  Real duration = args("duration", 1.0);
98 
99  for (unsigned int t_step = 0; t_step < duration/solid_system.deltat; t_step++)
100  {
101  // Progress in current phase [0..1]
102  Real progress = t_step * solid_system.deltat / duration;
103  systems.parameters.set<Real>("progress") = progress;
104  systems.parameters.set<unsigned int>("step") = t_step;
105 
106  // Update message
107 
108  out << "===== Time Step " << std::setw(4) << t_step;
109  out << " (" << std::fixed << std::setprecision(2) << std::setw(6) << (progress*100.) << "%)";
110  out << ", time = " << std::setw(7) << solid_system.time;
111  out << " =====" << std::endl;
112 
113  // Advance timestep in auxiliary system
114  aux_system.current_local_solution->close();
115  aux_system.old_local_solution->close();
116  *aux_system.older_local_solution = *aux_system.old_local_solution;
117  *aux_system.old_local_solution = *aux_system.current_local_solution;
118 
119  out << "Solving Solid" << std::endl;
120  solid_system.solve();
121  aux_system.reinit();
122 
123  // Carry out the adaptive mesh refinement/coarsening
124  out << "Doing a reinit of the equation systems" << std::endl;
125  systems.reinit();
126 
127  if (t_step % args("output/frequency", 1) == 0)
128  {
129  std::stringstream file_name;
130  file_name << args("results_directory", "./")
131  << "fem_"
132  << std::setw(6)
133  << std::setfill('0')
134  << t_step
135  << ".pvtu";
136 
137  io->write_equation_systems(file_name.str(), systems);
138  }
139  // Advance to the next timestep in a transient problem
140  out << "Advancing to next step" << std::endl;
141  solid_system.time_solver->advance_timestep();
142  }
143 }
144 
145 
146 
147 int main(int argc, char ** argv)
148 {
149  // Initialize libMesh and any dependent libraries
150  LibMeshInit init(argc, argv);
151 
152  // This example requires a linear solver package.
153  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
154  "--enable-petsc, --enable-trilinos, or --enable-eigen");
155 
156  // Skip this example if we do not meet certain requirements
157 #ifndef LIBMESH_HAVE_VTK
158  libmesh_example_requires(false, "--enable-vtk");
159 #endif
160 
161  // Trilinos gives us an inverted element on this one...
162  libmesh_example_requires(libMesh::default_solver_package() != TRILINOS_SOLVERS, "--enable-petsc");
163 
164  // Threaded assembly doesn't currently work with the moving mesh
165  // code.
166  // We'll skip this example for now.
167  if (libMesh::n_threads() > 1)
168  {
169  libMesh::out << "We skip fem_system_ex2 when using threads.\n"
170  << std::endl;
171  return 0;
172  }
173 
174  // read simulation parameters from file
175  GetPot args = GetPot("solid.in");
176 
177  // Create System and Mesh
178  int dim = args("mesh/generation/dimension", 3);
179  libmesh_example_requires(dim <= LIBMESH_DIM, "3D support");
180 
181  // Create a mesh distributed across the default MPI communicator.
182  Mesh mesh(init.comm(), dim);
183 
184  EquationSystems systems(mesh);
185 
186  // Create and set systems up
187  setup(systems, mesh, args);
188 
189  // run the systems
190  run_timestepping(systems, args);
191 
192  out << "Finished calculations" << std::endl;
193  return 0;
194 }
void run_timestepping(EquationSystems &systems, GetPot &args)
This is the EquationSystems class.
UniquePtr< NumericVector< Number > > old_local_solution
All the values I need to compute my contribution to the simulation at hand.
unsigned int n_threads()
Definition: libmesh_base.h:125
unsigned int dim
ElemType
Defines an enum for geometric element types.
MeshBase & mesh
This class provides a specific system class.
virtual void reinit() libmesh_override
Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be use...
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.
This class implements reading and writing meshes in the VTK format.
Definition: vtk_io.h:59
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
SolverPackage default_solver_package()
Definition: libmesh.C:995
GetPot args
Definition: solid_system.h:59
virtual void reinit()
Reinitialize all the systems.
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.
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
virtual void reinit() libmesh_override
Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be use...
UniquePtr< NumericVector< Number > > older_local_solution
All the values I need to compute my contribution to the simulation at hand.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
int main(int argc, char **argv)
T & set(const std::string &)
Definition: parameters.h:469
OStreamProxy out
void setup(EquationSystems &systems, Mesh &mesh, GetPot &args)
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
void save_initial_mesh()
Definition: solid_system.C:61
Parameters parameters
Data structure holding arbitrary parameters.
const MeshBase & get_mesh() const
const T_sys & get_system(const std::string &name) 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 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...