Functions
int	main (int argc, char **argv)
Function Documentation

◆ main()

int main	(	int	argc,
		char **	argv
	)
Definition at line 72 of file vector_fe_ex9.C.
 {
   // Initialize libMesh.
   LibMeshInit init(argc, argv);
 
   // This example requires PETSc
   libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS,
                            "--enable-petsc");
 
   // Parse the input file.
   GetPot infile("vector_fe_ex9.in");
 
   // But allow the command line to override it.
   infile.parse_command_line(argc, argv);
 
   // Read in parameters from the command line and the input file.
   const unsigned int dimension = 2;
   const unsigned int grid_size = infile("grid_size", 2);
   const bool mms = infile("mms", true);
   const Real nu = infile("nu", 1.);
   const bool cavity = infile("cavity", false);
 
   // Skip higher-dimensional examples on a lower-dimensional libMesh build.
   libmesh_example_requires(dimension <= LIBMESH_DIM, dimension << "D support");
 
   // Create a mesh, with dimension to be overridden later, distributed
   // across the default MPI communicator.
   Mesh mesh(init.comm());
 
   // Use the MeshTools::Generation mesh generator to create a uniform
   // grid on the cube [-1,1]^D. To accomodate first order side hierarchics, we must
   // use TRI6/7 elements
   const std::string elem_str = infile("element_type", std::string("TRI6"));
 
   libmesh_error_msg_if(elem_str != "TRI6" && elem_str != "TRI7",
                        "You selected "
                            << elem_str
                            << " but this example must be run with TRI6, TRI7, QUAD8, or QUAD9 in 2d"
                            << " or with TET14, or HEX27 in 3d.");
 
   if (mms && !cavity)
     MeshTools::Generation::build_square(
         mesh, grid_size, grid_size, 0., 2., -1, 1., Utility::string_to_enum<ElemType>(elem_str));
   else if (cavity)
     MeshTools::Generation::build_square(
         mesh, grid_size, grid_size, -1, 1, -1, 1., Utility::string_to_enum<ElemType>(elem_str));
   else
     MeshTools::Generation::build_square(mesh,
                                         5 * grid_size,
                                         grid_size,
                                         0.,
                                         10.,
                                         -1,
                                         1.,
                                         Utility::string_to_enum<ElemType>(elem_str));
   mesh.print_info();
 
   // Create an equation systems object.
   EquationSystems equation_systems(mesh);
 
   // Declare the system "Mixed" and its variables.
   auto & system = equation_systems.add_system<NonlinearImplicitSystem>("HDG");
 
   // Adds the velocity variables and their gradients
   system.add_variable("qu", FIRST, L2_LAGRANGE_VEC);
   system.add_variable("qv", FIRST, L2_LAGRANGE_VEC);
   system.add_variable("vel_x", FIRST, L2_LAGRANGE);
   system.add_variable("vel_y", FIRST, L2_LAGRANGE);
 
   // Add our Lagrange multiplier to the implicit system
   system.add_variable("lm_u", FIRST, SIDE_HIERARCHIC);
   system.add_variable("lm_v", FIRST, SIDE_HIERARCHIC);
   const auto p_num = system.add_variable("pressure", FIRST, L2_LAGRANGE);
   if (cavity)
     system.add_variable("global_lm", FIRST, SCALAR);
 
   const FEType vector_fe_type(FIRST, L2_LAGRANGE_VEC);
   const FEType scalar_fe_type(FIRST, L2_LAGRANGE);
   const FEType lm_fe_type(FIRST, SIDE_HIERARCHIC);
 
   StaticCondensation * sc = nullptr;
   if (system.has_static_condensation())
     {
       sc = &system.get_static_condensation();
       sc->dont_condense_vars({p_num});
     }
 
   HDGProblem hdg(nu, cavity);
   hdg.mesh = &mesh;
   hdg.system = &system;
   hdg.dof_map = &system.get_dof_map();
   hdg.vector_fe = FEVectorBase::build(dimension, vector_fe_type);
   hdg.scalar_fe = FEBase::build(dimension, scalar_fe_type);
   hdg.qrule = QBase::build(QGAUSS, dimension, FIFTH);
   hdg.qface = QBase::build(QGAUSS, dimension - 1, FIFTH);
   hdg.vector_fe_face = FEVectorBase::build(dimension, vector_fe_type);
   hdg.scalar_fe_face = FEBase::build(dimension, scalar_fe_type);
   hdg.lm_fe_face = FEBase::build(dimension, lm_fe_type);
   hdg.mms = mms;
   hdg.sc = sc;
 
   system.nonlinear_solver->residual_object = &hdg;
   system.nonlinear_solver->jacobian_object = &hdg;
 
   hdg.init();
 
   // Initialize the data structures for the equation system.
   equation_systems.init();
   equation_systems.print_info();
 
   // Solve the implicit system for the Lagrange multiplier
   system.solve();
 
   if (mms)
   {
     //
     // Now we will compute our solution approximation errors
     //
 
     equation_systems.parameters.set<const ExactSoln *>("vel_x_exact_sol") = &hdg.u_true_soln;
     equation_systems.parameters.set<const ExactSoln *>("vel_y_exact_sol") = &hdg.v_true_soln;
     equation_systems.parameters.set<const ExactSoln *>("pressure_exact_sol") = &hdg.p_true_soln;
     ExactSolution exact_sol(equation_systems);
     exact_sol.attach_exact_value(&compute_error);
 
     // Compute the error.
     exact_sol.compute_error("HDG", "vel_x");
     exact_sol.compute_error("HDG", "vel_y");
     exact_sol.compute_error("HDG", "pressure");
 
     // Print out the error values.
     libMesh::out << "L2 error for u is: " << exact_sol.l2_error("HDG", "vel_x") << std::endl;
     libMesh::out << "L2 error for v is: " << exact_sol.l2_error("HDG", "vel_y") << std::endl;
     libMesh::out << "L2 error for pressure is: " << exact_sol.l2_error("HDG", "pressure")
                  << std::endl;
   }
 
 #ifdef LIBMESH_HAVE_EXODUS_API
 
   // We write the file in the ExodusII format.
   ExodusII_IO(mesh).write_equation_systems("out.e", equation_systems);
 
 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
 
   // All done.
   return 0;
 }
Functions

Function Documentation

◆ main()
