reduced_basis_ex5.C File Reference

Go to the source code of this file.

Functions
void	scale_mesh_and_plot (EquationSystems &es, const RBParameters &mu, const std::string &filename)
void	compute_stresses (EquationSystems &es)
int	main (int argc, char **argv)

Function Documentation

compute_stresses()

void compute_stresses

(

EquationSystems &

es

)

Definition at line 404 of file reduced_basis_ex5.C.

References libMesh::TypeVector< T >::add_scaled(), libMesh::FEGenericBase< OutputType >::build(), libMesh::System::current_solution(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), elasticity_tensor(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::System::n_vars(), libMesh::Utility::pow(), libMesh::DenseMatrix< T >::resize(), libMesh::DenseMatrix< T >::scale(), libMesh::System::solution, libMesh::System::update(), libMesh::System::variable_number(), and libMesh::DofMap::variable_type().

Referenced by scale_mesh_and_plot().

{
  LOG_SCOPE("compute_stresses()", "main");

  const MeshBase & mesh = es.get_mesh();

  const unsigned int dim = mesh.mesh_dimension();
  libmesh_assert_less_equal(dim, 3);

  ElasticityRBConstruction & system = es.get_system<ElasticityRBConstruction>("RBElasticity");

  unsigned int displacement_vars[3];
  displacement_vars[0] = system.variable_number ("u");
  if (dim > 1)
    displacement_vars[1] = system.variable_number ("v");
  if (dim > 2)
    displacement_vars[2] = system.variable_number ("w");
  const unsigned int u_var = system.variable_number ("u");

  const DofMap & dof_map = system.get_dof_map();
  FEType fe_type = dof_map.variable_type(u_var);
  std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
  QGauss qrule (dim, fe_type.default_quadrature_order());
  fe->attach_quadrature_rule (&qrule);

  const std::vector<Real> & JxW = fe->get_JxW();
  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();

  // Also, get a reference to the ExplicitSystem
  ExplicitSystem & stress_system = es.get_system<ExplicitSystem>("StressSystem");
  const DofMap & stress_dof_map = stress_system.get_dof_map();
  unsigned int sigma_vars[3][3];
  sigma_vars[0][0] = stress_system.variable_number ("sigma_00");

  if (dim > 1)
    {
      sigma_vars[0][1] = stress_system.variable_number ("sigma_01");
      sigma_vars[1][0] = stress_system.variable_number ("sigma_10");
      sigma_vars[1][1] = stress_system.variable_number ("sigma_11");
    }

  if (dim > 2)
    {
      sigma_vars[0][2] = stress_system.variable_number ("sigma_02");
      sigma_vars[1][2] = stress_system.variable_number ("sigma_12");
      sigma_vars[2][0] = stress_system.variable_number ("sigma_20");
      sigma_vars[2][1] = stress_system.variable_number ("sigma_21");
      sigma_vars[2][2] = stress_system.variable_number ("sigma_22");
    }
  unsigned int vonMises_var = stress_system.variable_number ("vonMises");

  // Storage for the stress dof indices on each element
  std::vector<std::vector<dof_id_type>> dof_indices_var(system.n_vars());
  std::vector<dof_id_type> stress_dof_indices_var;

  // To store the stress tensor on each element
  DenseMatrix<Number> elem_sigma;

  for (const auto & elem : mesh.active_local_element_ptr_range())
    {
      // We'll only be computing stresses on the primary manifold of a
      // mesh, not, for instance, on spline nodes
      if (elem->dim() != dim)
        continue;

      for (unsigned int var=0; var<dim; var++)
        dof_map.dof_indices (elem, dof_indices_var[var], displacement_vars[var]);

      fe->reinit (elem);

      elem_sigma.resize(dim, dim);

      for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        for (unsigned int C_i=0; C_i<dim; C_i++)
          for (unsigned int C_j=0; C_j<dim; C_j++)
            for (unsigned int C_k=0; C_k<dim; C_k++)
              {
                const unsigned int n_var_dofs = dof_indices_var[C_k].size();

                // Get the gradient at this quadrature point
                Gradient displacement_gradient;
                for (unsigned int l=0; l<n_var_dofs; l++)
                  displacement_gradient.add_scaled(dphi[l][qp], system.current_solution(dof_indices_var[C_k][l]));

                for (unsigned int C_l=0; C_l<dim; C_l++)
                  elem_sigma(C_i,C_j) += JxW[qp] * (elasticity_tensor(C_i, C_j, C_k, C_l) * displacement_gradient(C_l));
              }

      // Get the average stresses by dividing by the element volume
      elem_sigma.scale(1./elem->volume());

      // load elem_sigma data into stress_system
      for (unsigned int i=0; i<dim; i++)
        for (unsigned int j=0; j<dim; j++)
          {
            stress_dof_map.dof_indices (elem, stress_dof_indices_var, sigma_vars[i][j]);

            // We are using CONSTANT MONOMIAL basis functions, hence we only need to get
            // one dof index per variable
            dof_id_type dof_index = stress_dof_indices_var[0];

            if ((stress_system.solution->first_local_index() <= dof_index) &&
                (dof_index < stress_system.solution->last_local_index()))
              {
                stress_system.solution->set(dof_index, elem_sigma(i,j));
              }

          }

      // Also, the von Mises stress

      // We're solving with general plane stress if in 2D, or with
      // uniaxial stress if in 1D
      Number sigma00 = elem_sigma(0,0);
      Number sigma11 = (dim > 1) ? elem_sigma(1,1) : 0;
      Number sigma22 = (dim > 2) ? elem_sigma(2,2) : 0;
      Number sum_term = pow(sigma00 - sigma11,2);
      if (dim > 1)
        sum_term += pow(sigma11 - sigma22,2) +
                    6.*pow(elem_sigma(0,1),2);
      if (dim > 2)
        sum_term += pow(sigma22 - sigma00,2) +
                    6.*(pow(elem_sigma(1,2),2) + pow(elem_sigma(2,0),2));
      Number vonMises_value = std::sqrt(0.5*sum_term);
      stress_dof_map.dof_indices (elem, stress_dof_indices_var, vonMises_var);
      dof_id_type dof_index = stress_dof_indices_var[0];
      if ((stress_system.solution->first_local_index() <= dof_index) &&
          (dof_index < stress_system.solution->last_local_index()))
        {
          stress_system.solution->set(dof_index, vonMises_value);
        }

    }

  // Should call close and update when we set vector entries directly
  stress_system.solution->close();
  stress_system.update();
}

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 86 of file reduced_basis_ex5.C.

References libMesh::BoundaryInfo::add_node(), libMesh::BoundaryInfo::add_side(), libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), libMesh::MeshTools::Generation::build_cube(), libMesh::CG, libMesh::ParallelObject::comm(), libMesh::command_line_next(), libMesh::CONSTANT, libMesh::default_solver_package(), dim, libMesh::MeshBase::get_boundary_info(), libMesh::LinearImplicitSystem::get_linear_solver(), libMesh::RBConstruction::get_rb_evaluation(), libMesh::BoundaryInfo::has_boundary_id(), libMesh::HEX27, libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::RBConstruction::initialize_rb_construction(), libMesh::INVALID_SOLVER_PACKAGE, libMesh::RBEvaluation::legacy_write_offline_data_to_files(), libMesh::RBConstruction::load_rb_solution(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::MONOMIAL, libMesh::RBConstruction::print_info(), libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::RBConstruction::process_parameters_file(), libMesh::MeshBase::read(), libMesh::RBDataDeserialization::RBEvaluationDeserialization::read_from_file(), libMesh::Real, scale_mesh_and_plot(), libMesh::RBConstruction::set_rb_evaluation(), libMesh::LinearSolver< T >::set_solver_type(), libMesh::RBParameters::set_value(), libMesh::RBConstruction::train_reduced_basis(), libMesh::TRILINOS_SOLVERS, ElasticityRBConstruction::var_order, libMesh::RBEvaluation::write_out_basis_functions(), and libMesh::RBDataSerialization::RBEvaluationSerialization::write_to_file().

{
  // Initialize libMesh.
  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");

#if !defined(LIBMESH_HAVE_XDR)
  // We need XDR support to write out reduced bases
  libmesh_example_requires(false, "--enable-xdr");
#elif defined(LIBMESH_DEFAULT_SINGLE_PRECISION)
  // XDR binary support requires double precision
  libmesh_example_requires(false, "--disable-singleprecision");
#elif defined(LIBMESH_DEFAULT_TRIPLE_PRECISION)
  // I have no idea why long double isn't working here... [RHS]
  libmesh_example_requires(false, "double precision");
#endif

  // Trilinos reports "true residual is too large" in the offline mode
  // portion of this example
  libmesh_example_requires(libMesh::default_solver_package() != TRILINOS_SOLVERS, "--enable-petsc or --enable-laspack");

  // This example only works on all meshes if libMesh was compiled for 3D
  libmesh_example_requires(LIBMESH_DIM == 3, "3D support");

#ifndef LIBMESH_ENABLE_DIRICHLET
  libmesh_example_requires(false, "--enable-dirichlet");
#else

  std::string parameters_filename = "reduced_basis_ex5.in";
  GetPot infile(parameters_filename);

  // Override input file arguments from the command line
  infile.parse_command_line(argc, argv);

  unsigned int n_elem_x = infile("n_elem_x", 0);
  unsigned int n_elem_y = infile("n_elem_y", 0);
  unsigned int n_elem_z = infile("n_elem_z", 0);
  Real x_size           = infile("x_size", 0.);
  Real y_size           = infile("y_size", 0.);
  Real z_size           = infile("z_size", 0.);

  bool store_basis_functions = infile("store_basis_functions", true);

  // Read the "online_mode" flag from the command line
  const int online_mode = libMesh::command_line_next("-online_mode", 0);

  Mesh mesh (init.comm());

  // Read a mesh from file?
  const std::string mesh_filename = infile("mesh_filename", "NO MESH FILE");

  if (mesh_filename == "NO MESH FILE")
    MeshTools::Generation::build_cube
      (mesh, n_elem_x, n_elem_y, n_elem_z,
       0., x_size, 0., y_size, 0., z_size, HEX27);
  else
    {
      mesh.read(mesh_filename);

      // We might be using legacy reduced_basis I/O, which relies on
      // Hilbert curve renumbering ... and if we have overlapping
      // nodes, as in the case of IGA meshes where FE nodes and spline
      // nodes can coincide, then we need unique_id() values to
      // disambiguate them when sorting.
#if !defined(LIBMESH_HAVE_CAPNPROTO) && !defined(LIBMESH_ENABLE_UNIQUE_ID)
  libmesh_example_requires(false, "--enable-unique-id or --enable-capnp");
#endif

      // We don't yet support the experimental BEXT sideset spec, so
      // for now we'll add our required sidesets manually for our BEXT
      // file.

      // To avoid extraneous allocation when obtaining element side vertex averages
      ElemSideBuilder side_builder;

      // Each processor should know about each boundary condition it can
      // see, so we loop over all elements, not just local elements.
      for (const auto & elem : mesh.element_ptr_range())
        for (auto side : elem->side_index_range())
          if (!elem->neighbor_ptr(side))
            {
              const Point side_center = side_builder(*elem, side).vertex_average();
              // Yes, BOUNDARY_ID_MIN_X is a weird ID to use at
              // min(z), but we got an IGA cantilever mesh with the
              // lever arm in the z direction
              if (side_center(2) < 0.1)
                mesh.get_boundary_info().add_side(elem, side, BOUNDARY_ID_MIN_X);
              if (side_center(2) > 11.9)
                mesh.get_boundary_info().add_side(elem, side, BOUNDARY_ID_MAX_X);
            }
    }

  // Let's add a node boundary condition so that we can impose a point
  // load and get more test coverage in the build_cube case.

  // Each processor should know about each boundary condition it can
  // see, so we loop over all elements, not just local elements.
  for (const auto & elem : mesh.element_ptr_range())
    {
      unsigned int side_max_x = 0, side_max_y = 0, side_max_z = 0;
      bool found_side_max_x = false, found_side_max_y = false, found_side_max_z = false;
      for (auto side : elem->side_index_range())
        {
          if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_X))
            {
              side_max_x = side;
              found_side_max_x = true;
            }

          if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_Y))
            {
              side_max_y = side;
              found_side_max_y = true;
            }

          if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_Z))
            {
              side_max_z = side;
              found_side_max_z = true;
            }
        }

      // If elem has sides on boundaries
      // BOUNDARY_ID_MAX_X, BOUNDARY_ID_MAX_Y, BOUNDARY_ID_MAX_Z
      // then let's set a node boundary condition
      if (found_side_max_x && found_side_max_y && found_side_max_z)
        {
          for (auto n : elem->node_index_range())
            {
              if (elem->is_node_on_side(n, side_max_x) &&
                  elem->is_node_on_side(n, side_max_y) &&
                  elem->is_node_on_side(n, side_max_z))
                {
                  mesh.get_boundary_info().add_node(elem->node_ptr(n), NODE_BOUNDARY_ID);
                }
            }
        }
    }

  mesh.print_info();

  const unsigned int dim = mesh.mesh_dimension();

  // Create an equation systems object.
  EquationSystems equation_systems(mesh);

  // We override RBConstruction with ElasticityRBConstruction in order to
  // specialize a few functions for this particular problem.
  ElasticityRBConstruction & rb_con =
    equation_systems.add_system<ElasticityRBConstruction>("RBElasticity");

  unsigned int order = infile("order", 1);

  rb_con.var_order = Order(order);

  // Also, initialize an ExplicitSystem to store stresses
  ExplicitSystem & stress_system =
    equation_systems.add_system<ExplicitSystem> ("StressSystem");
  stress_system.add_variable("sigma_00", CONSTANT, MONOMIAL);

  // Lots of ifs to keep order consistent
  if (dim > 1)
    stress_system.add_variable("sigma_01", CONSTANT, MONOMIAL);
  if (dim > 2)
      stress_system.add_variable("sigma_02", CONSTANT, MONOMIAL);
  if (dim > 1)
      stress_system.add_variable("sigma_10", CONSTANT, MONOMIAL);
  if (dim > 1)
      stress_system.add_variable("sigma_11", CONSTANT, MONOMIAL);
  if (dim > 2)
    {
      stress_system.add_variable("sigma_12", CONSTANT, MONOMIAL);
      stress_system.add_variable("sigma_20", CONSTANT, MONOMIAL);
      stress_system.add_variable("sigma_21", CONSTANT, MONOMIAL);
      stress_system.add_variable("sigma_22", CONSTANT, MONOMIAL);
    }
  stress_system.add_variable("vonMises", CONSTANT, MONOMIAL);

  // Initialize the data structures for the equation system.
  equation_systems.init ();
  equation_systems.print_info();

  // Build a new RBEvaluation object which will be used to perform
  // Reduced Basis calculations. This is required in both the
  // "Offline" and "Online" stages.
  ElasticityRBEvaluation rb_eval(mesh.comm());

  // We need to give the RBConstruction object a pointer to
  // our RBEvaluation object
  rb_con.set_rb_evaluation(rb_eval);

  if (!online_mode) // Perform the Offline stage of the RB method
    {
      // Use CG solver.  This echoes the Petsc command line options set
      // in run.sh, but also works for non-PETSc linear solvers.
      rb_con.get_linear_solver()->set_solver_type(CG);

      // Read in the data that defines this problem from the specified text file
      rb_con.process_parameters_file(parameters_filename);

      // Print out info that describes the current setup of rb_con
      rb_con.print_info();

      // Prepare rb_con for the Construction stage of the RB method.
      // This sets up the necessary data structures and performs
      // initial assembly of the "truth" affine expansion of the PDE.
      rb_con.initialize_rb_construction();

      // Compute the reduced basis space by computing "snapshots", i.e.
      // "truth" solves, at well-chosen parameter values and employing
      // these snapshots as basis functions.
      rb_con.train_reduced_basis();

      // Write out the data that will subsequently be required for the Evaluation stage
#if defined(LIBMESH_HAVE_CAPNPROTO)
      RBDataSerialization::RBEvaluationSerialization rb_eval_writer(rb_con.get_rb_evaluation());
      rb_eval_writer.write_to_file("rb_eval.bin");
#else
      rb_con.get_rb_evaluation().legacy_write_offline_data_to_files();
#endif

      // If requested, write out the RB basis functions for visualization purposes
      if (store_basis_functions)
        {
          // Write out the basis functions
          rb_con.get_rb_evaluation().write_out_basis_functions(rb_con);
        }
    }
  else // Perform the Online stage of the RB method
    {
      // Read in the reduced basis data
#if defined(LIBMESH_HAVE_CAPNPROTO)
      RBDataDeserialization::RBEvaluationDeserialization rb_eval_reader(rb_eval);
      rb_eval_reader.read_from_file("rb_eval.bin", /*read_error_bound_data*/ true);
#else
      rb_eval.legacy_read_offline_data_from_files();
#endif

      // Initialize online parameters
      Real online_x_scaling = infile("online_x_scaling", 0.);
      Real online_load_Fx   = infile("online_load_Fx",   0.);
      Real online_load_Fy   = infile("online_load_Fy",   0.);
      Real online_load_Fz   = infile("online_load_Fz",   0.);
      Real online_point_load_Fx   = infile("online_point_load_Fx",   0.);
      Real online_point_load_Fy   = infile("online_point_load_Fy",   0.);
      Real online_point_load_Fz   = infile("online_point_load_Fz",   0.);
      RBParameters online_mu;
      online_mu.set_value("x_scaling", online_x_scaling);
      online_mu.set_value("load_Fx",   online_load_Fx);
      online_mu.set_value("load_Fy",   online_load_Fy);
      online_mu.set_value("load_Fz",   online_load_Fz);
      online_mu.set_value("point_load_Fx",   online_point_load_Fx);
      online_mu.set_value("point_load_Fy",   online_point_load_Fy);
      online_mu.set_value("point_load_Fz",   online_point_load_Fz);
      rb_eval.set_parameters(online_mu);
      rb_eval.print_parameters();

      // Now do the Online solve using the precomputed reduced basis
      rb_eval.rb_solve(rb_eval.get_n_basis_functions());

      if (store_basis_functions)
        {
          // Read in the basis functions
          rb_eval.read_in_basis_functions(rb_con);

          // Plot the solution
          rb_con.load_rb_solution();

          const RBParameters & rb_eval_params = rb_eval.get_parameters();
          scale_mesh_and_plot(equation_systems, rb_eval_params, "RB_sol");
        }
    }

#endif // LIBMESH_ENABLE_DIRICHLET

  return 0;
}

scale_mesh_and_plot()

void scale_mesh_and_plot	(	EquationSystems &	es,
		const RBParameters &	mu,
		const std::string &	filename
	)

Definition at line 369 of file reduced_basis_ex5.C.

References libMesh::ParallelObject::comm(), compute_stresses(), libMesh::MeshBase::get_constraint_rows(), libMesh::EquationSystems::get_mesh(), libMesh::RBParameters::get_value(), libMesh::MeshBase::is_serial(), libMesh::libmesh_ignore(), mesh, TIMPI::Communicator::min(), libMesh::MeshOutput< MT >::write_equation_systems(), and libMesh::ExodusII_IO::write_equation_systems().

Referenced by main().

{
  // Loop over the mesh nodes and move them!
  MeshBase & mesh = es.get_mesh();

  for (auto & node : mesh.node_ptr_range())
    (*node)(0) *= mu.get_value("x_scaling");

  // Post-process the solution to compute the stresses
  compute_stresses(es);

#ifdef LIBMESH_HAVE_EXODUS_API
  // Distributed IGA meshes don't yet support re-gathering to serial
  bool do_exodus_write =
    (es.get_mesh().get_constraint_rows().empty() ||
     es.get_mesh().is_serial());
  es.comm().min(do_exodus_write);
  if (do_exodus_write)
    ExodusII_IO (mesh).write_equation_systems (file_basename+".e", es);
#ifdef LIBMESH_HAVE_NEMESIS_API
  else
    Nemesis_IO (mesh).write_equation_systems (file_basename+".nem", es);
#endif
#endif

  // Avoid unused variable warnings when Exodus is not enabled
  libmesh_ignore(file_basename);

  // Loop over the mesh nodes and move them!
  for (auto & node : mesh.node_ptr_range())
    (*node)(0) /= mu.get_value("x_scaling");
}