fem_system_ex2.C File Reference

Go to the source code of this file.

Functions
void	setup (EquationSystems &systems, Mesh &mesh, GetPot &args)
void	run_timestepping (EquationSystems &systems, GetPot &args)
int	main (int argc, char **argv)

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 157 of file fem_system_ex2.C.

References libMesh::default_solver_package(), dim, libMesh::TriangleWrapper::init(), libMesh::INVALID_SOLVER_PACKAGE, mesh, libMesh::n_threads(), libMesh::out, run_timestepping(), setup(), and libMesh::TRILINOS_SOLVERS.

{
  // Initialize libMesh and any dependent libraries
  LibMeshInit init(argc, argv);

  // This example requires a linear solver package.
  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
                           "--enable-petsc, --enable-trilinos, or --enable-eigen");

  // Trilinos gives us an inverted element on this one...
  libmesh_example_requires(libMesh::default_solver_package() != TRILINOS_SOLVERS, "--enable-petsc");

  // Threaded assembly doesn't currently work with the moving mesh
  // code.
  // We'll skip this example for now.
  if (libMesh::n_threads() > 1)
    {
      libMesh::out << "We skip fem_system_ex2 when using threads.\n"
                   << std::endl;
      return 0;
    }

  // read simulation parameters from file
  GetPot args = GetPot("solid.in");

  // But allow the command line to override
  args.parse_command_line(argc, argv);

  // Create System and Mesh
  int dim = args("mesh/generation/dimension", 3);
  libmesh_example_requires(dim <= LIBMESH_DIM, "3D support");

  // Create a mesh distributed across the default MPI communicator.
  Mesh mesh(init.comm(), cast_int<unsigned char>(dim));

  EquationSystems systems(mesh);

  // Create and set systems up
  setup(systems, mesh, args);

  // run the systems
  run_timestepping(systems, args);

  out << "Finished calculations" << std::endl;
  return 0;
}

run_timestepping()

void run_timestepping	(	EquationSystems &	systems,
		GetPot &	args
	)

Definition at line 93 of file fem_system_ex2.C.

References libMesh::NumericVector< T >::close(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::TransientSystem< Base >::old_local_solution, libMesh::TransientSystem< Base >::older_local_solution, libMesh::out, libMesh::EquationSystems::parameters, libMesh::Real, libMesh::EquationSystems::reinit(), and libMesh::Parameters::set().

Referenced by main().

{
  TransientExplicitSystem & aux_system = systems.get_system<TransientExplicitSystem>("auxiliary");

  SolidSystem & solid_system = systems.get_system<SolidSystem>("solid");

#ifdef LIBMESH_HAVE_VTK
  std::unique_ptr<VTKIO> io = std::make_unique<VTKIO>(systems.get_mesh());
#endif

  Real duration = args("duration", 1.0);

  for (unsigned int t_step = 0; t_step < duration/solid_system.deltat; t_step++)
    {
      // Progress in current phase [0..1]
      Real progress = t_step * solid_system.deltat / duration;
      systems.parameters.set<Real>("progress") = progress;
      systems.parameters.set<unsigned int>("step") = t_step;

      // Update message

      out << "===== Time Step " << std::setw(4) << t_step;
      out << " (" << std::fixed << std::setprecision(2) << std::setw(6) << (progress*100.) << "%)";
      out << ", time = " << std::setw(7) << solid_system.time;
      out << " =====" << std::endl;

      // Advance timestep in auxiliary system
      aux_system.current_local_solution->close();
      aux_system.old_local_solution->close();
      *aux_system.older_local_solution = *aux_system.old_local_solution;
      *aux_system.old_local_solution = *aux_system.current_local_solution;

      out << "Solving Solid" << std::endl;
      solid_system.solve();
      aux_system.current_local_solution->close();
      *aux_system.solution = *aux_system.current_local_solution;
      aux_system.solution->close();

      // Carry out the adaptive mesh refinement/coarsening
      out << "Doing a reinit of the equation systems" << std::endl;
      systems.reinit();

#ifdef LIBMESH_HAVE_VTK
      if (t_step % args("output/frequency", 1) == 0)
        {
          std::stringstream file_name;
          file_name << args("results_directory", "./")
                    << "fem_"
                    << std::setw(6)
                    << std::setfill('0')
                    << t_step
                    << ".pvtu";

          io->write_equation_systems(file_name.str(), systems);
        }
#endif
      // Advance to the next timestep in a transient problem
      out << "Advancing to next step" << std::endl;
      solid_system.time_solver->advance_timestep();
    }
}

setup()

void setup	(	EquationSystems &	systems,
		Mesh &	mesh,
		GetPot &	args
	)

Definition at line 49 of file fem_system_ex2.C.

References libMesh::EquationSystems::add_system(), SolidSystem::args, libMesh::MeshTools::Generation::build_cube(), libMesh::System::current_local_solution, dim, libMesh::EquationSystems::init(), libMesh::libmesh_assert(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::EquationSystems::parameters, libMesh::MeshBase::print_info(), SolidSystem::save_initial_mesh(), libMesh::Parameters::set(), and libMesh::System::solution.

Referenced by main().

{
  const unsigned int dim = mesh.mesh_dimension();
  // We currently invert tensors with the assumption that they're 3x3
  libmesh_assert (dim == 3);

  // Generating Mesh
  ElemType eltype = Utility::string_to_enum<ElemType>(args("mesh/generation/element_type", "hex8"));
  int nx = args("mesh/generation/num_elem", 4, 0);
  int ny = args("mesh/generation/num_elem", 4, 1);
  int nz = dim > 2 ? args("mesh/generation/num_elem", 4, 2) : 0;
  double origx = args("mesh/generation/origin", -1.0, 0);
  double origy = args("mesh/generation/origin", -1.0, 1);
  double origz = args("mesh/generation/origin", 0.0, 2);
  double sizex = args("mesh/generation/size", 2.0, 0);
  double sizey = args("mesh/generation/size", 2.0, 1);
  double sizez = args("mesh/generation/size", 2.0, 2);
  MeshTools::Generation::build_cube(mesh, nx, ny, nz,
                                    origx, origx+sizex, origy, origy+sizey, origz, origz+sizez, eltype);
  mesh.print_info();

  // Creating Systems
  SolidSystem & imms = systems.add_system<SolidSystem> ("solid");
  imms.args = args;

  // Build up auxiliary system
  ExplicitSystem & aux_sys = systems.add_system<TransientExplicitSystem>("auxiliary");

  // Initialize the system
  systems.parameters.set<unsigned int>("phase") = 0;
  systems.init();

  imms.save_initial_mesh();

  // Fill global solution vector from local ones
  aux_sys.current_local_solution->close();
  *aux_sys.solution = *aux_sys.current_local_solution;
  aux_sys.solution->close();
}