vector_fe_ex5.C File Reference

Go to the source code of this file.

Namespaces
namespace	libMesh
	The libMesh namespace provides an interface to certain functionality in the library.

Functions
void	compute_residual (const NumericVector< Number > &X, NumericVector< Number > &R, NonlinearImplicitSystem &S)
void	compute_jacobian (const NumericVector< Number > &X, SparseMatrix< Number > &J, NonlinearImplicitSystem &S)
int	main (int argc, char **argv)

Function Documentation

compute_jacobian()

void compute_jacobian	(	const NumericVector< Number > &	X,
		SparseMatrix< Number > &	J,
		NonlinearImplicitSystem &	S
	)

Definition at line 315 of file assembly.C.

{
  // It is a good idea to make sure we are assembling
  // the proper system.
  libmesh_assert_equal_to(system.name(), "Poisson");
 
  // Get the DG parameters
  auto & es = system.get_equation_systems();
  const auto sigma = es.parameters.get<Real>("sigma");
  const auto epsilon = es.parameters.get<Real>("epsilon");
 
  // Get a constant reference to the mesh object.
  const MeshBase & mesh = system.get_mesh();
 
  // The dimension that we are running
  const auto dim = mesh.mesh_dimension();
 
  // A reference to the  DofMap object for this system.  The  DofMap
  // object handles the index translation from node and element numbers
  // to degree of freedom numbers.  We will talk more about the  DofMap
  // in future examples.
  const DofMap & dof_map = system.get_dof_map();
 
  // Get a constant reference to the Finite Element type
  // for the first (and only) variable in the system.
  FEType fe_type = dof_map.variable_type(0);
 
  // Build a Finite Element object of the specified type.
  // Note that FEVectorBase is a typedef for the templated FE
  // class.
  std::unique_ptr<FEVectorBase> fe(FEVectorBase::build(dim, fe_type));
 
  // An automatically determined Gauss quadrature rule for numerical integration.
  QGauss qrule(dim, fe_type.default_quadrature_order());
 
  // Tell the finite element object to use our quadrature rule.
  fe->attach_quadrature_rule(&qrule);
 
  // Declare a special finite element object for face integration.
  std::unique_ptr<FEVectorBase> fe_face(FEVectorBase::build(dim, fe_type));
  // And for neighbor integration
  std::unique_ptr<FEVectorBase> fe_neighbor_face(FEVectorBase::build(dim, fe_type));
 
  // Boundary integration requires one quadrature rule,
  // with dimensionality one less than the dimensionality
  // of the element.
  QGauss qface(dim - 1, fe_type.default_quadrature_order());
 
  // Tell the face finite element objects to use our
  // quadrature rule.
  fe_face->attach_quadrature_rule(&qface);
  fe_neighbor_face->attach_quadrature_rule(&qface);
 
  // Here we define some references to cell-specific data that
  // will be used to assemble the linear system.
  //
  // The element Jacobian * quadrature weight at each integration point.
  const auto & JxW = fe->get_JxW();
 
  // The element shape function gradients evaluated at the quadrature points.
  const auto & dphi = fe->get_dphi();
 
  // face integration points
  const auto & face_xyz = fe_face->get_xyz();
 
  // Face shape function values
  const auto & phi_face = fe_face->get_phi();
 
  // Face shape function gradients
  const auto & dphi_face = fe_face->get_dphi();
 
  // Neighbor shape function values
  const auto & phi_neighbor = fe_neighbor_face->get_phi();
 
  // Neighbor face shape function gradients
  const auto & dphi_neighbor = fe_neighbor_face->get_dphi();
 
  // face normals
  const auto & normals = fe_face->get_normals();
 
  // face JxW
  const auto & JxW_face = fe_face->get_JxW();
 
  DenseMatrix<Number> Kee;
  DenseMatrix<Number> Ken;
  DenseMatrix<Number> Kne;
  DenseMatrix<Number> Knn;
 
  // This vector will hold the degree of freedom indices for
  // the element.  These define where in the global system
  // the element degrees of freedom get mapped.
  std::vector<dof_id_type> dof_indices;
  std::vector<dof_id_type> dof_indices_neighbor;
 
  // To avoid extraneous allocation when building element sides (to compute their volumes)
  ElemSideBuilder side_builder;
 
  // Now we will loop over all the elements in the mesh.
  // We will compute the element matrix and right-hand-side
  // contribution.
  //
  // Element iterators are a nice way to iterate through all the
  // elements, or all the elements that have some property.  The
  // iterator el will iterate from the first to the last element on
  // the local processor.  The iterator end_el tells us when to stop.
  // It is smart to make this one const so that we don't accidentally
  // mess it up!  In case users later modify this program to include
  // refinement, we will be safe and will only consider the active
  // elements; hence we use a variant of the active_elem_iterator.
  for (const auto & elem : mesh.active_local_element_ptr_range())
  {
    // Get the degree of freedom indices for the
    // current element.  These define where in the global
    // matrix and right-hand-side this element will
    // contribute to.
    dof_map.dof_indices(elem, dof_indices);
 
    // Compute the element-specific data for the current
    // element.  This involves computing the location of the
    // quadrature points (q_point) and the shape functions
    // (phi, dphi) for the current element.
    fe->reinit(elem);
 
    libmesh_assert_msg(dphi.size() == dof_indices.size(),
                       "dphi size doesn't match dof_indices size");
 
    // DenseMatrix::resize() member will automatically zero out the matrix.
    Kee.resize(dof_indices.size(), dof_indices.size());
 
    // diffusion elemental jacobian
    for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
      for (std::size_t i = 0; i < dof_indices.size(); i++)
        for (std::size_t j = 0; j < dof_indices.size(); j++)
          Kee(i, j) += JxW[qp] * dphi[i][qp].contract(dphi[j][qp]);
 
    // Now we consider jacobian contributions from the sides
    for (auto side : elem->side_index_range())
    {
      // We need to compute h for penalty terms
      const auto side_volume = side_builder(*elem, side).volume();
      fe_face->reinit(elem, side);
      const auto elem_b_order = static_cast<unsigned int>(fe_face->get_order());
      const auto h_elem = elem->volume() / side_volume / std::pow(elem_b_order, 2);
 
      // No neighbor means we must be on a boundary
      if (!elem->neighbor_ptr(side))
      {
        for (unsigned int qp = 0; qp < qface.n_points(); qp++)
          for (std::size_t i = 0; i < dof_indices.size(); i++)
            for (std::size_t j = 0; j < dof_indices.size(); j++)
 
            {
              Kee(i, j) -= dphi_face[j][qp] * normals[qp] * phi_face[i][qp] * JxW_face[qp];
              Kee(i, j) +=
                  epsilon * phi_face[j][qp] * dphi_face[i][qp] * normals[qp] * JxW_face[qp];
              Kee(i, j) += sigma / h_elem * phi_face[j][qp] * phi_face[i][qp] * JxW_face[qp];
            }
      }
      else // We must be on an interior side
      {
        const Elem * neighbor = elem->neighbor_ptr(side);
 
        const auto elem_id = elem->id();
        const auto neighbor_id = neighbor->id();
 
        // We don't want to erroneously add multiple contributions from the same interior face
        if ((neighbor->active() && (neighbor->level() == elem->level()) &&
             (elem_id < neighbor_id)) ||
            (neighbor->level() < elem->level()))
        {
          dof_map.dof_indices(neighbor, dof_indices_neighbor);
 
          // Make sure we have the matching quadrature points on face and neighbor
          std::vector<Point> neighbor_xyz;
          FEMap::inverse_map(elem->dim(), neighbor, face_xyz,
                             neighbor_xyz);
 
          fe_neighbor_face->reinit(neighbor, &neighbor_xyz);
 
          libmesh_assert_msg(dphi_neighbor.size() == dof_indices_neighbor.size(),
                             "dphi_neighbor size doesn't match dof_indices_neighbor size");
 
          Ken.resize(dof_indices.size(), dof_indices_neighbor.size());
          Kne.resize(dof_indices_neighbor.size(), dof_indices.size());
          Knn.resize(dof_indices_neighbor.size(), dof_indices_neighbor.size());
 
          // Now add the DG contribution to the local jacobian
          for (unsigned int qp = 0; qp < qface.n_points(); qp++)
          {
            // element-element contribution
            for (std::size_t i = 0; i < dof_indices.size(); i++)
              for (std::size_t j = 0; j < dof_indices.size(); j++)
              {
                Kee(i, j) -= 0.5 * dphi_face[j][qp] * normals[qp] * phi_face[i][qp] * JxW_face[qp];
                Kee(i, j) +=
                    epsilon * 0.5 * phi_face[j][qp] * dphi_face[i][qp] * normals[qp] * JxW_face[qp];
                Kee(i, j) += sigma / h_elem * phi_face[j][qp] * phi_face[i][qp] * JxW_face[qp];
              }
            // element-neighbor contribution
            for (std::size_t i = 0; i < dof_indices.size(); i++)
              for (std::size_t j = 0; j < dof_indices_neighbor.size(); j++)
              {
                Ken(i, j) -=
                    0.5 * dphi_neighbor[j][qp] * normals[qp] * phi_face[i][qp] * JxW_face[qp];
                Ken(i, j) += epsilon * 0.5 * -phi_neighbor[j][qp] * dphi_face[i][qp] * normals[qp] *
                             JxW_face[qp];
                Ken(i, j) += sigma / h_elem * -phi_neighbor[j][qp] * phi_face[i][qp] * JxW_face[qp];
              }
            // Neighbor-element contribution
            for (std::size_t i = 0; i < dof_indices_neighbor.size(); i++)
              for (std::size_t j = 0; j < dof_indices_neighbor.size(); j++)
              {
                Kne(i, j) +=
                    0.5 * dphi_face[j][qp] * normals[qp] * phi_neighbor[i][qp] * JxW_face[qp];
                Kne(i, j) -= epsilon * 0.5 * phi_face[j][qp] * dphi_neighbor[i][qp] * normals[qp] *
                             JxW_face[qp];
                Kne(i, j) -= sigma / h_elem * phi_face[j][qp] * phi_neighbor[i][qp] * JxW_face[qp];
              }
            // Neighbor-neighbor contribution
            for (std::size_t i = 0; i < dof_indices_neighbor.size(); i++)
              for (std::size_t j = 0; j < dof_indices_neighbor.size(); j++)
              {
                Knn(i, j) +=
                    0.5 * dphi_neighbor[j][qp] * normals[qp] * phi_neighbor[i][qp] * JxW_face[qp];
                Knn(i, j) -= epsilon * 0.5 * -phi_neighbor[j][qp] * dphi_neighbor[i][qp] *
                             normals[qp] * JxW_face[qp];
                Knn(i, j) -=
                    sigma / h_elem * -phi_neighbor[j][qp] * phi_neighbor[i][qp] * JxW_face[qp];
              }
          }
 
          J.add_matrix(Ken, dof_indices, dof_indices_neighbor);
          J.add_matrix(Kne, dof_indices_neighbor, dof_indices);
          J.add_matrix(Knn, dof_indices_neighbor, dof_indices_neighbor);
        } // whether we've done this internal side integral before
      }   // whether we're on an internal side
    }     // loop over sides
 
    J.add_matrix(Kee, dof_indices, dof_indices);
  } // element loop
}

References libMesh::Elem::active(), libMesh::SparseMatrix< T >::add_matrix(), libMesh::FEGenericBase< OutputType >::build(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::System::get_mesh(), libMesh::DofObject::id(), libMesh::FEMap::inverse_map(), libMesh::Elem::level(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::QBase::n_points(), libMesh::System::name(), libMesh::Elem::neighbor_ptr(), libMesh::EquationSystems::parameters, libMesh::Real, libMesh::DenseMatrix< T >::resize(), and libMesh::DofMap::variable_type().

Referenced by libMesh::FirstOrderUnsteadySolver::compute_second_order_eqns(), LaplaceSystem::element_time_derivative(), PoissonSystem::element_time_derivative(), HeatSystem::element_time_derivative(), SigmaPhysics::element_time_derivative(), main(), and LaplaceSystem::side_constraint().

compute_residual()

void compute_residual	(	const NumericVector< Number > &	X,
		NumericVector< Number > &	R,
		NonlinearImplicitSystem &	S
	)

Vectors to hold the local solution degree of freedom values

Vector to hold the local solution

Definition at line 27 of file assembly.C.

{
  // It is a good idea to make sure we are assembling
  // the proper system.
  libmesh_assert_equal_to(system.name(), "Poisson");
 
  // Get the DG parameters
  auto & es = system.get_equation_systems();
  const auto sigma = es.parameters.get<Real>("sigma");
  const auto epsilon = es.parameters.get<Real>("epsilon");
 
  // Get a constant reference to the mesh object.
  const MeshBase & mesh = system.get_mesh();
 
  // The dimension that we are running
  const auto dim = mesh.mesh_dimension();
 
  // A reference to the  DofMap object for this system.  The  DofMap
  // object handles the index translation from node and element numbers
  // to degree of freedom numbers.  We will talk more about the  DofMap
  // in future examples.
  const DofMap & dof_map = system.get_dof_map();
 
  // Get a constant reference to the Finite Element type
  // for the first (and only) variable in the system.
  FEType fe_type = dof_map.variable_type(0);
 
  // Build a Finite Element object of the specified type.
  // Note that FEVectorBase is a typedef for the templated FE
  // class.
  std::unique_ptr<FEVectorBase> fe(FEVectorBase::build(dim, fe_type));
 
  // An automatically determined Gauss quadrature rule for numerical integration.
  QGauss qrule(dim, fe_type.default_quadrature_order());
 
  // Tell the finite element object to use our quadrature rule.
  fe->attach_quadrature_rule(&qrule);
 
  // Declare a special finite element object for face integration.
  std::unique_ptr<FEVectorBase> fe_face(FEVectorBase::build(dim, fe_type));
  // And for neighbor integration
  std::unique_ptr<FEVectorBase> fe_neighbor_face(FEVectorBase::build(dim, fe_type));
 
  // Boundary integration requires one quadrature rule,
  // with dimensionality one less than the dimensionality
  // of the element.
  QGauss qface(dim - 1, fe_type.default_quadrature_order());
 
  // Tell the face finite element objects to use our
  // quadrature rule.
  fe_face->attach_quadrature_rule(&qface);
  fe_neighbor_face->attach_quadrature_rule(&qface);
 
  // Here we define some references to cell-specific data that
  // will be used to assemble the linear system.
  //
  // The element Jacobian * quadrature weight at each integration point.
  const auto & JxW = fe->get_JxW();
 
  // element integration points
  const auto & xyz = fe->get_xyz();
 
  // The element shape function values evaluated at the quadrature points
  const auto & phi = fe->get_phi();
 
  // The element shape function gradients evaluated at the quadrature points.
  const auto & dphi = fe->get_dphi();
 
  // face integration points
  const auto & face_xyz = fe_face->get_xyz();
 
  // Face shape function values
  const auto & phi_face = fe_face->get_phi();
 
  // Face shape function gradients
  const auto & dphi_face = fe_face->get_dphi();
 
  // Neighbor shape function values
  const auto & phi_neighbor = fe_neighbor_face->get_phi();
 
  // Neighbor face shape function gradients
  const auto & dphi_neighbor = fe_neighbor_face->get_dphi();
 
  // face normals
  const auto & normals = fe_face->get_normals();
 
  // face JxW
  const auto & JxW_face = fe_face->get_JxW();
 
  DenseVector<Number> Fe;
  DenseVector<Number> Fn;
 
  // This vector will hold the degree of freedom indices for
  // the element.  These define where in the global system
  // the element degrees of freedom get mapped.
  std::vector<dof_id_type> dof_indices;
  std::vector<dof_id_type> dof_indices_neighbor;
 
  std::vector<Number> dof_u;
  std::vector<Number> dof_u_neighbor;
 
  std::vector<VectorValue<Number>> u;
  std::vector<TensorValue<Number>> grad_u;
  std::vector<VectorValue<Number>> u_neighbor;
  std::vector<TensorValue<Number>> grad_u_neighbor;
 
  // Now we will loop over all the elements in the mesh.
  // We will compute the element matrix and right-hand-side
  // contribution.
  //
  // Element iterators are a nice way to iterate through all the
  // elements, or all the elements that have some property.  The
  // iterator el will iterate from the first to the last element on
  // the local processor.  The iterator end_el tells us when to stop.
  // It is smart to make this one const so that we don't accidentally
  // mess it up!  In case users later modify this program to include
  // refinement, we will be safe and will only consider the active
  // elements; hence we use a variant of the active_elem_iterator.
  for (const auto & elem : mesh.active_local_element_ptr_range())
  {
    // Get the degree of freedom indices for the
    // current element.  These define where in the global
    // matrix and right-hand-side this element will
    // contribute to.
    dof_map.dof_indices(elem, dof_indices);
 
    // Compute the element-specific data for the current
    // element.  This involves computing the location of the
    // quadrature points (q_point) and the shape functions
    // (phi, dphi) for the current element.
    fe->reinit(elem);
 
    libmesh_assert_msg(dphi.size() == dof_indices.size(),
                       "dphi size doesn't match dof_indices size");
 
    // DenseVector::resize() member will automatically zero out the vector.
    Fe.resize(dof_indices.size());
 
    // Get the local solution vector
    X.get(dof_indices, dof_u);
 
    // build the element solution gradient
    grad_u.resize(qrule.n_points());
    for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
    {
      grad_u[qp] = 0;
      for (std::size_t i = 0; i < dof_indices.size(); i++)
        grad_u[qp] += dof_u[i] * dphi[i][qp];
    }
 
    for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
    {
      auto forcing_func = VectorValue<Number>(forcing_function(0, xyz[qp](0), xyz[qp](1)),
                                              forcing_function(1, xyz[qp](0), xyz[qp](1)),
                                              0);
      for (std::size_t i = 0; i < dof_indices.size(); i++)
      {
        // diffusion elemental residual
        Fe(i) += JxW[qp] * dphi[i][qp].contract(grad_u[qp]);
 
        // forcing function
        Fe(i) += JxW[qp] * phi[i][qp] * forcing_func;
      }
    }
 
    // To avoid extraneous allocation when building element sides (to compute their volumes)
    ElemSideBuilder side_builder;
 
    // Now we consider residual contributions from the sides
    for (auto side : elem->side_index_range())
    {
      // We need to compute h for penalty terms
      const auto side_volume = side_builder(*elem, side).volume();
      fe_face->reinit(elem, side);
      const auto elem_b_order = static_cast<unsigned int>(fe_face->get_order());
      const auto h_elem = elem->volume() / side_volume / std::pow(elem_b_order, 2);
 
      // build the face solution value and gradient
      u.resize(qface.n_points());
      grad_u.resize(qface.n_points());
      for (unsigned int qp = 0; qp < qface.n_points(); qp++)
      {
        u[qp] = 0;
        grad_u[qp] = 0;
        for (std::size_t i = 0; i < dof_indices.size(); i++)
        {
          u[qp] += dof_u[i] * phi_face[i][qp];
          grad_u[qp] += dof_u[i] * dphi_face[i][qp];
        }
      }
 
      // No neighbor means we must be on a boundary
      if (!elem->neighbor_ptr(side))
      {
        for (unsigned int qp = 0; qp < qface.n_points(); qp++)
        {
          auto fn = VectorValue<Number>(exact_solution(0, face_xyz[qp](0), face_xyz[qp](1)),
                                        exact_solution(1, face_xyz[qp](0), face_xyz[qp](1)),
                                        0);
          for (std::size_t i = 0; i < dof_indices.size(); i++)
          {
            Fe(i) -= grad_u[qp] * normals[qp] * phi_face[i][qp] * JxW_face[qp];
            Fe(i) += epsilon * (u[qp] - fn) * dphi_face[i][qp] * normals[qp] * JxW_face[qp];
            Fe(i) += sigma / h_elem * (u[qp] - fn) * phi_face[i][qp] * JxW_face[qp];
          }
        }
      }
      else // We must be on an interior side
      {
        const Elem * neighbor = elem->neighbor_ptr(side);
 
        const auto elem_id = elem->id();
        const auto neighbor_id = neighbor->id();
 
        // We don't want to erroneously add multiple contributions from the same interior face
        if ((neighbor->active() && (neighbor->level() == elem->level()) &&
             (elem_id < neighbor_id)) ||
            (neighbor->level() < elem->level()))
        {
          dof_map.dof_indices(neighbor, dof_indices_neighbor);
 
          // Make sure we have the matching quadrature points on face and neighbor
          std::vector<Point> neighbor_xyz;
          FEMap::inverse_map(elem->dim(), neighbor, face_xyz,
                             neighbor_xyz);
 
          fe_neighbor_face->reinit(neighbor, &neighbor_xyz);
 
          libmesh_assert_msg(dphi_neighbor.size() == dof_indices_neighbor.size(),
                             "dphi_neighbor size doesn't match dof_indices_neighbor size");
 
          Fn.resize(dof_indices_neighbor.size());
          X.get(dof_indices_neighbor, dof_u_neighbor);
 
          // build the neighbor solution value and gradient
          u_neighbor.resize(qface.n_points());
          grad_u_neighbor.resize(qface.n_points());
          for (unsigned int qp = 0; qp < qface.n_points(); qp++)
          {
            u_neighbor[qp] = 0;
            grad_u_neighbor[qp] = 0;
            for (std::size_t i = 0; i < dof_indices_neighbor.size(); i++)
            {
              u_neighbor[qp] += dof_u_neighbor[i] * phi_neighbor[i][qp];
              grad_u_neighbor[qp] += dof_u_neighbor[i] * dphi_neighbor[i][qp];
            }
          }
 
          // Now add the DG contribution to the local residual
          for (unsigned int qp = 0; qp < qface.n_points(); qp++)
          {
            // element contribution
            for (std::size_t i = 0; i < dof_indices.size(); i++)
            {
              Fe(i) -= 0.5 * (grad_u[qp] * normals[qp] + grad_u_neighbor[qp] * normals[qp]) *
                       phi_face[i][qp] * JxW_face[qp];
              Fe(i) += epsilon * 0.5 * (u[qp] - u_neighbor[qp]) * dphi_face[i][qp] * normals[qp] *
                       JxW_face[qp];
              Fe(i) += sigma / h_elem * (u[qp] - u_neighbor[qp]) * phi_face[i][qp] * JxW_face[qp];
            }
            // Neighbor contribution
            for (std::size_t i = 0; i < dof_indices_neighbor.size(); i++)
            {
              Fn(i) += 0.5 * (grad_u[qp] * normals[qp] + grad_u_neighbor[qp] * normals[qp]) *
                       phi_neighbor[i][qp] * JxW_face[qp];
              Fn(i) -= epsilon * 0.5 * (u[qp] - u_neighbor[qp]) * dphi_neighbor[i][qp] *
                       normals[qp] * JxW_face[qp];
              Fn(i) -=
                  sigma / h_elem * (u[qp] - u_neighbor[qp]) * phi_neighbor[i][qp] * JxW_face[qp];
            }
          }
 
          R.add_vector(Fn, dof_indices_neighbor);
        } // whether we've done this internal side integral before
      }   // whether we're on an internal side
    }     // loop over sides
 
    R.add_vector(Fe, dof_indices);
  } // element loop
}

References libMesh::Elem::active(), libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), exact_solution, forcing_function, libMesh::NumericVector< T >::get(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::System::get_mesh(), libMesh::DofObject::id(), libMesh::FEMap::inverse_map(), libMesh::Elem::level(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::QBase::n_points(), libMesh::System::name(), libMesh::Elem::neighbor_ptr(), libMesh::EquationSystems::parameters, libMesh::Real, libMesh::DenseVector< T >::resize(), and libMesh::DofMap::variable_type().

Referenced by main().

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 63 of file vector_fe_ex5.C.

{
  // Initialize libraries.
  LibMeshInit init(argc, argv);
 
  // This example requires a non-linear solver package.
  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS ||
                               libMesh::default_solver_package() == TRILINOS_SOLVERS,
                           "--enable-petsc or --enable-trilinos");
 
  // Parse input file
  GetPot input_file("vector_fe_ex5.in");
 
  // But allow the command line to override it.
  input_file.parse_command_line(argc, argv);
 
  // Read DG parameters
  const Real epsilon = input_file("epsilon", -1);
  const Real sigma = input_file("sigma", 6);
 
  // Read mesh size
  const std::size_t nx = input_file("nx", 15);
  const std::size_t ny = input_file("ny", 15);
 
  // Read approximation order - default to FIRST
  const unsigned int approx_order = input_file("approx_order", 1);
 
  // Brief message to the user regarding the program name
  // and command line arguments.
  libMesh::out << "Running " << argv[0];
 
  for (int i = 1; i < argc; i++)
    libMesh::out << " " << argv[i];
 
  libMesh::out << std::endl << std::endl;
 
  // Skip this 2D example if libMesh was compiled as 1D-only.
  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
 
  // Create a mesh, with dimension to be overridden later, on the
  // default MPI communicator.
  Mesh mesh(init.comm());
 
  // Use the MeshTools::Generation mesh generator to create a uniform
  // 2D grid on the square [-1,1]^2.  We instruct the mesh generator
  // to build a mesh of 15x15 QUAD4 elements. Note that at the end of this
  // function call, the mesh will be prepared and remote elements will be deleted
  MeshTools::Generation::build_square(mesh, nx, ny, -1., 1., -1., 1., QUAD4);
 
  // Print information about the mesh to the screen.
  mesh.print_info();
 
  // Create an equation systems object.
  EquationSystems equation_systems(mesh);
  equation_systems.parameters.set<Real>("epsilon") = epsilon;
  equation_systems.parameters.set<Real>("sigma") = sigma;
 
  // Declare the Poisson system and its variables.
  // The Poisson system is another example of a steady system.
  auto & nl_system = equation_systems.add_system<NonlinearImplicitSystem>("Poisson");
 
  // Adds the variable "u" to "Poisson". Variable "u" will be approximated using
  // first-order (can be overridden with command line) vector monomial elements.
  // Since the mesh is 2-D, "u" will have two components.
  nl_system.add_variable("u", static_cast<Order>(approx_order), MONOMIAL_VEC);
 
  // Set the residual and Jacobian evaluation functions
  auto & nl_solver = *nl_system.nonlinear_solver;
  nl_solver.residual = compute_residual;
  nl_solver.jacobian = compute_jacobian;
 
  // Initialize the data structures for the equation system.
  equation_systems.init();
 
  // Prints information about the system to the screen.
  equation_systems.print_info();
 
  // A p=0 "solve" here is a hack to test ExodusII output; the
  // solve() is stymied by division by zero.
  std::unique_ptr<FPEDisabler> maybe_disable_fpe;
  if (! approx_order) /* CONSTANT, MONOMIAL_VEC */
    maybe_disable_fpe = std::make_unique<FPEDisabler>();
 
  nl_system.solve();
 
#ifdef LIBMESH_HAVE_EXODUS_API
  if (! approx_order) /* CONSTANT, MONOMIAL_VEC */
    {
      // Warn about failed solve for CONSTANT approximations (Poisson needs at least p=1)
      libMesh::out << "\nWARNING: The elemental vector order is CONSTANT. The solution" << std::endl
                   << "written out to 'out_constant.e' will be trivial." << std::endl;
 
      // Need to get string of variable names (vector components) for 'set_output_variables()'
      std::vector<std::string> varnames;
      equation_systems.build_variable_names(varnames);
 
      // Create ExodusII object with an unambiguous filename (indicating this solution is special)
      ExodusII_IO exo_ptr(mesh);
      exo_ptr.write("out_constant.e");
      exo_ptr.set_output_variables(varnames);
      exo_ptr.write_element_data(equation_systems);
    }
  else
    ExodusII_IO(mesh).write_equation_systems("out.e", equation_systems);
#endif
 
  // All done.
  return 0;
}

References libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), libMesh::MeshTools::Generation::build_square(), libMesh::EquationSystems::build_variable_names(), compute_jacobian(), compute_residual(), libMesh::default_solver_package(), libMesh::EquationSystems::init(), libMesh::SparseMatrix< T >::init(), main(), mesh, libMesh::MONOMIAL_VEC, libMesh::out, libMesh::EquationSystems::parameters, libMesh::PETSC_SOLVERS, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::QUAD4, libMesh::Real, libMesh::Parameters::set(), libMesh::ExodusII_IO::set_output_variables(), libMesh::TRILINOS_SOLVERS, libMesh::ExodusII_IO::write(), libMesh::ExodusII_IO::write_element_data(), and libMesh::ExodusII_IO::write_equation_systems().