miscellaneous_ex1.C File Reference

Go to the source code of this file.

Functions
void	assemble_wave (EquationSystems &es, const std::string &system_name)
int	main (int argc, char **argv)
void	assemble_wave (EquationSystems &es, const std::string &libmesh_dbg_var(system_name))

Function Documentation

assemble_wave() [1/2]

void assemble_wave	(	EquationSystems &	es,
		const std::string &	system_name
	)

Definition at line 318 of file transient_ex2.C.

References libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), dim, libMesh::DofMap::dof_indices(), libMesh::Parameters::get(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::QBase::n_points(), libMesh::EquationSystems::parameters, libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), and libMesh::SECOND.

Referenced by main().

{
  // It is a good idea to make sure we are assembling
  // the proper system.
  libmesh_assert_equal_to (system_name, "Wave");

  // Get a constant reference to the mesh object.
  const MeshBase & mesh = es.get_mesh();

  // The dimension that we are running.
  const unsigned int dim = mesh.mesh_dimension();

  // Copy the speed of sound to a local variable.
  const Real speed = es.parameters.get<Real>("speed");

  // If we added Neumann conditions we would need density too
  // const Real rho   = es.parameters.get<Real>("fluid density");

  // Get a reference to our system, as before.
  NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);

  // Get a constant reference to the Finite Element type
  // for the first (and only) variable in the system.
  FEType fe_type = t_system.get_dof_map().variable_type(0);

  // In here, we will add the element matrices to the
  // @e additional matrices "stiffness_mass" and "damping"
  // and the additional vector "force", not to the members
  // "matrix" and "rhs".  Therefore, get writable
  // references to them.
  SparseMatrix<Number> & stiffness = t_system.get_matrix("stiffness");
  SparseMatrix<Number> & damping   = t_system.get_matrix("damping");
  SparseMatrix<Number> & mass      = t_system.get_matrix("mass");
  NumericVector<Number> & force    = t_system.get_vector("force");

  // Some solver packages (PETSc) are especially picky about
  // allocating sparsity structure and truly assigning values
  // to this structure.  Namely, matrix additions, as performed
  // later, exhibit acceptable performance only for identical
  // sparsity structures.  Therefore, explicitly zero the
  // values in the collective matrix, so that matrix additions
  // encounter identical sparsity structures.
  SparseMatrix<Number> & matrix = *t_system.matrix;
  DenseMatrix<Number> zero_matrix;

  // Build a Finite Element object of the specified type.  Since the
  // FEBase::build() member dynamically creates memory we will
  // store the object as a std::unique_ptr<FEBase>.  This can be thought
  // of as a pointer that will clean up after itself.
  std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));

  // A 2nd order Gauss quadrature rule for numerical integration.
  QGauss qrule (dim, SECOND);

  // Tell the finite element object to use our quadrature rule.
  fe->attach_quadrature_rule (&qrule);

  // The element Jacobian * quadrature weight at each integration point.
  const std::vector<Real> & JxW = fe->get_JxW();

  // The element shape functions evaluated at the quadrature points.
  const std::vector<std::vector<Real>> & phi = fe->get_phi();

  // The element shape function gradients evaluated at the quadrature
  // points.
  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();

  // 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.
  const DofMap & dof_map = t_system.get_dof_map();

  // The element mass, damping and stiffness matrices
  // and the element contribution to the rhs.
  DenseMatrix<Number> Ke, Ce, Me;
  DenseVector<Number> Fe;

  // 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;

  // Now we will loop over all the elements in the mesh.
  // We will compute the element matrix and right-hand-side
  // contribution.
  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);

      // Zero the element matrices and rhs before
      // summing them.  We use the resize member here because
      // the number of degrees of freedom might have changed from
      // the last element.  Note that this will be the case if the
      // element type is different (i.e. the last element was HEX8
      // and now have a PRISM6).
      {
        const unsigned int n_dof_indices = dof_indices.size();

        Ke.resize          (n_dof_indices, n_dof_indices);
        Ce.resize          (n_dof_indices, n_dof_indices);
        Me.resize          (n_dof_indices, n_dof_indices);
        zero_matrix.resize (n_dof_indices, n_dof_indices);
        Fe.resize          (n_dof_indices);
      }

      // Now loop over the quadrature points.  This handles
      // the numeric integration.
      for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        {
          // Now we will build the element matrix.  This involves
          // a double loop to integrate the test functions (i) against
          // the trial functions (j).
          for (std::size_t i=0; i<phi.size(); i++)
            for (std::size_t j=0; j<phi.size(); j++)
              {
                Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
                Me(i,j) += JxW[qp]*phi[i][qp]*phi[j][qp]
                  *1./(speed*speed);
              } // end of the matrix summation loop
        } // end of quadrature point loop

      // Now compute the contribution to the element matrix and the
      // right-hand-side vector if the current element lies on the
      // boundary.
      {
        // In this example no natural boundary conditions will
        // be considered.  The code is left here so it can easily
        // be extended.
        //
        // don't do this for any side
#if 0
        for (auto side : elem->side_index_range())
          if (elem->neighbor_ptr(side) == nullptr)
            {
              // Declare a special finite element object for
              // boundary integration.
              std::unique_ptr<FEBase> fe_face (FEBase::build(dim, fe_type));

              // Boundary integration requires one quadrature rule,
              // with dimensionality one less than the dimensionality
              // of the element.
              QGauss qface(dim-1, SECOND);

              // Tell the finite element object to use our
              // quadrature rule.
              fe_face->attach_quadrature_rule (&qface);

              // The value of the shape functions at the quadrature
              // points.
              const std::vector<std::vector<Real>> &  phi_face = fe_face->get_phi();

              // The Jacobian * Quadrature Weight at the quadrature
              // points on the face.
              const std::vector<Real> & JxW_face = fe_face->get_JxW();

              // Compute the shape function values on the element
              // face.
              fe_face->reinit(elem, side);

              // Here we consider a normal acceleration acc_n=1 applied to
              // the whole boundary of our mesh.
              const Real acc_n_value = 1.0;

              // Loop over the face quadrature points for integration.
              for (unsigned int qp=0; qp<qface.n_points(); qp++)
                {
                  // Right-hand-side contribution due to prescribed
                  // normal acceleration.
                  for (std::size_t i=0; i<phi_face.size(); i++)
                    {
                      Fe(i) += acc_n_value*rho
                        *phi_face[i][qp]*JxW_face[qp];
                    }
                } // end face quadrature point loop
            } // end if (elem->neighbor_ptr(side) == nullptr)
#endif // 0

        // In this example the Dirichlet boundary conditions will be
        // imposed via penalty method after the
        // system is assembled.

      } // end boundary condition section

      // If this assembly program were to be used on an adaptive mesh,
      // we would have to apply any hanging node constraint equations
      // by uncommenting the following lines:
      // std::vector<unsigned int> dof_indicesC = dof_indices;
      // std::vector<unsigned int> dof_indicesM = dof_indices;
      // dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
      // dof_map.constrain_element_matrix (Ce, dof_indicesC);
      // dof_map.constrain_element_matrix (Me, dof_indicesM);

      // Finally, simply add the contributions to the additional
      // matrices and vector.
      stiffness.add_matrix (Ke, dof_indices);
      damping.add_matrix   (Ce, dof_indices);
      mass.add_matrix      (Me, dof_indices);

      force.add_vector     (Fe, dof_indices);

      // For the overall matrix, explicitly zero the entries where
      // we added values in the other ones, so that we have
      // identical sparsity footprints.
      matrix.add_matrix(zero_matrix, dof_indices);

    } // end of element loop
}

assemble_wave() [2/2]

void assemble_wave	(	EquationSystems &	es,
		const std::string &	libmesh_dbg_varsystem_name
	)

Definition at line 239 of file miscellaneous_ex1.C.

References libMesh::DenseMatrix< T >::add(), libMesh::NumericVector< T >::add(), libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), libMesh::FEGenericBase< OutputType >::build_InfFE(), dim, libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::ImplicitSystem::get_system_matrix(), libMesh::libmesh_ignore(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::FEInterface::n_dofs(), libMesh::EquationSystems::parameters, libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::ExplicitSystem::rhs, libMesh::SECOND, libMesh::TOLERANCE, and libMesh::MeshTools::weight().

{
  // It is a good idea to make sure we are assembling
  // the proper system.
  libmesh_assert_equal_to (system_name, "Wave");

  // Avoid unused variable warnings when compiling without infinite
  // elements enabled.
  libmesh_ignore(es);

#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS

  // Get a constant reference to the mesh object.
  const MeshBase & mesh = es.get_mesh();

  // Get a reference to the system we are solving.
  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Wave");

  // 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.
  const DofMap & dof_map = system.get_dof_map();

  // The dimension that we are running.
  const unsigned int dim = mesh.mesh_dimension();

  // Copy the speed of sound to a local variable.
  const Real speed = es.parameters.get<Real>("speed");

  // Get a constant reference to the Finite Element type
  // for the first (and only) variable in the system.
  const FEType & fe_type = dof_map.variable_type(0);

  // Build a Finite Element object of the specified type.  Since the
  // FEBase::build() member dynamically creates memory we will
  // store the object as a std::unique_ptr<FEBase>.
  std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));

  // Do the same for an infinite element.
  std::unique_ptr<FEBase> inf_fe (FEBase::build_InfFE(dim, fe_type));

  // A 2nd order Gauss quadrature rule for numerical integration.
  QGauss qrule (dim, SECOND);

  // Tell the finite element object to use our quadrature rule.
  fe->attach_quadrature_rule (&qrule);

  // Due to its internal structure, the infinite element handles
  // quadrature rules differently.  It takes the quadrature
  // rule which has been initialized for the FE object, but
  // creates suitable quadrature rules by @e itself.  The user
  // need not worry about this.
  inf_fe->attach_quadrature_rule (&qrule);

  // Define data structures to contain the element matrix
  // and right-hand-side vector contribution.  Following
  // basic finite element terminology we will denote these
  // "Ke",  "Ce", "Me", and "Fe" for the stiffness, damping
  // and mass matrices, and the load vector.  Note that in
  // Acoustics, these descriptors though do @e not match the
  // true physical meaning of the projectors.  The final
  // overall system, however, resembles the conventional
  // notation again.
  DenseMatrix<Number> Ke;
  DenseMatrix<Number> Ce;
  DenseMatrix<Number> Me;
  DenseVector<Number> Fe;

  // 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;

  // The global system matrix
  SparseMatrix<Number> & matrix = system.get_system_matrix();

  // Now we will loop over all the elements in the mesh.
  // We will compute the element matrix and right-hand-side
  // contribution.
  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);

      const unsigned int n_dofs =
        cast_int<unsigned int>(dof_indices.size());

      // The mesh contains both finite and infinite elements.  These
      // elements are handled through different classes, namely
      // FE and InfFE, respectively.  However, since both
      // are derived from FEBase, they share the same interface,
      // and overall burden of coding is @e greatly reduced through
      // using a pointer, which is adjusted appropriately to the
      // current element type.
      FEBase * cfe = nullptr;

      // This here is almost the only place where we need to
      // distinguish between finite and infinite elements.
      // For faster computation, however, different approaches
      // may be feasible.
      //
      // Up to now, we do not know what kind of element we
      // have.  Aske the element of what type it is:
      if (elem->infinite())
        {
          // We have an infinite element.  Let cfe point
          // to our InfFE object.  This is handled through
          // a std::unique_ptr.  Through the std::unique_ptr::get() we "borrow"
          // the pointer, while the  std::unique_ptr inf_fe is
          // still in charge of memory management.
          cfe = inf_fe.get();
        }
      else
        {
          // This is a conventional finite element.  Let fe handle it.
          cfe = fe.get();

          // Boundary conditions.
          // Here we just zero the rhs-vector. For natural boundary
          // conditions check e.g. previous examples.
          {
            // Zero the RHS for this element.
            Fe.resize (n_dofs);

            system.rhs->add_vector (Fe, dof_indices);
          } // end boundary condition section
        } // else if (elem->infinite())

      // This is slightly different from the Poisson solver:
      // Since the finite element object may change, we have to
      // initialize the constant references to the data fields
      // each time again, when a new element is processed.
      //
      // Due to infinite extent of the element, the 'pure' Jacobian is divergent
      // and we must use a re-scaled one. As re-scaling, a 'decay'-function
      // is used (in 3D this is \f$ r^{-2} \f$).
      //
      // To account for this extra weight, \p phi, \p dphi and \p weight are
      // re-scaled as well:
      //  * J      --> J x decay^2
      //  * phi    --> phi/decay x r
      //  * dphi   --> dphi/decay x r
      //  * weight --> weight / r^2
      // With this, the product of Jacobian, test and trial functions (and their derivatives)
      // can be used as normal since the extra weights cancel.
      //
      // The element Jacobian * quadrature weight at each integration point.
      const std::vector<Real> & JxW = cfe->get_JxWxdecay_sq();

      // The element shape functions evaluated at the quadrature points.
      const std::vector<std::vector<Real>> & phi = cfe->get_phi_over_decayxR();
      //const std::vector<std::vector<Real>> & phi = cfe->get_phi();

      // The element shape function gradients evaluated at the quadrature
      // points, divided by weight x rad
      const std::vector<std::vector<RealGradient>> & dphi = cfe->get_dphi_over_decayxR();
      //const std::vector<std::vector<RealGradient>> & dphi = cfe->get_dphi();

      // The infinite elements need more data fields than conventional FE.
      // These are the gradients of the phase term dphase, an additional
      // radial weight for the test functions Sobolev_weight, and its
      // gradient.
      //
      // Note that these data fields are also initialized appropriately by
      // the FE method, so that the weak form (below) is valid for @e both
      // finite and infinite elements.
      const std::vector<RealGradient> & dphase  = cfe->get_dphase();
      const std::vector<Real> &         weight  = cfe->get_Sobolev_weightxR_sq();
      const std::vector<RealGradient> & dweight = cfe->get_Sobolev_dweightxR_sq();

      // Now this is all independent of whether we use an FE
      // or an InfFE.  Nice, hm? ;-)
      //
      // Compute the element-specific data, as described
      // in previous examples.
      cfe->reinit (elem);

      // Zero the element matrices.  Boundary conditions were already
      // processed in the FE-only section, see above.
      Ke.resize (n_dofs, n_dofs);
      Ce.resize (n_dofs, n_dofs);
      Me.resize (n_dofs, n_dofs);

      // The total number of quadrature points for infinite elements
      // @e has to be determined in a different way, compared to
      // conventional finite elements.  This type of access is also
      // valid for finite elements, so this can safely be used
      // anytime, instead of asking the quadrature rule, as
      // seen in previous examples.
      unsigned int max_qp = cfe->n_quadrature_points();

      // Loop over the quadrature points.
      for (unsigned int qp=0; qp<max_qp; qp++)
        {
          // Similar to the modified access to the number of quadrature
          // points, the number of shape functions may also be obtained
          // in a different manner.  This offers the great advantage
          // of being valid for both finite and infinite elements.
          const unsigned int n_sf =
            FEInterface::n_dofs(cfe->get_fe_type(), elem);

          // Now we will build the element matrices.  Since the infinite
          // elements are based on a Petrov-Galerkin scheme, the
          // resulting system matrices are non-symmetric. The additional
          // weight, described before, is part of the trial space.
          //
          // For the finite elements, though, these matrices are symmetric
          // just as we know them, since the additional fields dphase,
          // weight, and dweight are initialized appropriately.
          //
          // test functions:    weight[qp]*phi[i][qp]
          // trial functions:   phi[j][qp]
          // phase term:        phase[qp]
          //
          // derivatives are similar, but note that these are of type
          // Point, not of type Real.
          for (unsigned int i=0; i<n_sf; i++)
            for (unsigned int j=0; j<n_sf; j++)
              {
                // (ndt*Ht + nHt*d) * nH
                Ke(i,j) +=
                  (
                   (dweight[qp] * phi[i][qp] // Point * Real  = Point
                    +                        // +
                    dphi[i][qp] * weight[qp] // Point * Real  = Point
                    ) * dphi[j][qp]
                   ) * JxW[qp];

                // (d*Ht*nmut*nH - ndt*nmu*Ht*H - d*nHt*nmu*H)
                Ce(i,j) +=
                  (
                   (dphase[qp] * dphi[j][qp]) // (Point * Point) = Real
                   * weight[qp] * phi[i][qp]  // * Real * Real   = Real
                   -                          // -
                   (dweight[qp] * dphase[qp]) // (Point * Point) = Real
                   * phi[i][qp] * phi[j][qp]  // * Real * Real   = Real
                   -                          // -
                   (dphi[i][qp] * dphase[qp]) // (Point * Point) = Real
                   * weight[qp] * phi[j][qp]  // * Real * Real   = Real
                   ) * JxW[qp];

                // (d*Ht*H * (1 - nmut*nmu))
                Me(i,j) +=
                  (
                   (1. - (dphase[qp] * dphase[qp]))       // (Real  - (Point * Point)) = Real
                   * phi[i][qp] * phi[j][qp] * weight[qp] // * Real *  Real  * Real    = Real
                   ) * JxW[qp];

              } // end of the matrix summation loop
        } // end of quadrature point loop

      // The element matrices are now built for this element.
      // Collect them in Ke, and then add them to the global matrix.
      // The SparseMatrix::add_matrix() member does this for us.
      Ke.add(1./speed        , Ce);
      Ke.add(1./(speed*speed), Me);

      // If this assembly program were to be used on an adaptive mesh,
      // we would have to apply any hanging node constraint equations
      dof_map.constrain_element_matrix(Ke, dof_indices);

      matrix.add_matrix (Ke, dof_indices);
    } // end of element loop

  // Note that we have not applied any boundary conditions so far.
  // Here we apply a unit load at the node located at (0,0,0).
  for (const auto & node : mesh.local_node_ptr_range())
    if (std::abs((*node)(0)) < TOLERANCE &&
        std::abs((*node)(1)) < TOLERANCE &&
        std::abs((*node)(2)) < TOLERANCE)
      {
        // The global number of the respective degree of freedom.
        unsigned int dn = node->dof_number(0,0,0);

        system.rhs->add (dn, 1.);
      }

#else

  // dummy assert
  libmesh_assert_not_equal_to (es.get_mesh().mesh_dimension(), 1);

#endif //ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
}

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 86 of file miscellaneous_ex1.C.

References libMesh::EquationSystems::add_system(), assemble_wave(), libMesh::MeshTools::Generation::build_cube(), libMesh::InfElemBuilder::build_inf_elem(), libMesh::MeshBase::find_neighbors(), libMesh::FIRST, libMesh::EquationSystems::get_system(), libMesh::HEX8, libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), mesh, libMesh::out, libMesh::EquationSystems::parameters, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::Real, libMesh::EquationSystems::reinit(), libMesh::Parameters::set(), libMesh::MeshRefinement::uniformly_refine(), libMesh::WRITE, libMesh::ExodusII_IO::write(), and libMesh::EquationSystems::write().

{
  // Initialize libMesh, like in example 2.
  LibMeshInit init (argc, argv);

  // This example requires Infinite Elements
#ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
  libmesh_example_requires(false, "--enable-ifem");
#else

  // Skip this 3D example if libMesh was compiled as 1D/2D-only.
  libmesh_example_requires(3 <= LIBMESH_DIM, "3D support");

  // Create a serialized mesh, distributed across the default MPI
  // communicator.
  Mesh mesh(init.comm());

  // Get command line arguments for mesh size
  GetPot input(argc, argv);

  const unsigned int nx = input("nx", 4),
                     ny = input("ny", 4),
                     nz = input("nz", 4);

  // Use the internal mesh generator to create elements
  // on the square [-1,1]^3, of type Hex8.
  MeshTools::Generation::build_cube (mesh,
                                     nx, ny, nz,
                                     -1., 1.,
                                     -1., 1.,
                                     -1., 1.,
                                     HEX8);

  // Print information about the mesh to the screen.
  mesh.print_info();

  // Write the mesh before the infinite elements are added
#ifdef LIBMESH_HAVE_EXODUS_API
  ExodusII_IO(mesh).write ("orig_mesh.e");
#endif

  // Normally, when a mesh is imported or created in
  // libMesh, only conventional elements exist.  The infinite
  // elements used here, however, require prescribed
  // nodal locations (with specified distances from an imaginary
  // origin) and configurations that a conventional mesh creator
  // in general does not offer.  Therefore, an efficient method
  // for building infinite elements is offered.  It can account
  // for symmetry planes and creates infinite elements in a fully
  // automatic way.
  //
  // Right now, the simplified interface is used, automatically
  // determining the origin.  Check MeshBase for a generalized
  // method that can even return the element faces of interior
  // vibrating surfaces.  The bool determines whether to be
  // verbose.
  InfElemBuilder builder(mesh);
  builder.build_inf_elem(true);

  // Reassign subdomain_id() of all infinite elements.
  // Otherwise, the exodus-api will fail on the mesh.
  for (auto & elem : mesh.element_ptr_range())
    if (elem->infinite())
      elem->subdomain_id() = 1;

  // Print information about the mesh to the screen.
  mesh.print_info();

  // Write the mesh with the infinite elements added.
  // Compare this to the original mesh.
#ifdef LIBMESH_HAVE_EXODUS_API
  ExodusII_IO(mesh).write ("ifems_added.e");
#endif

  // After building infinite elements, we have to let
  // the elements find their neighbors again.
  mesh.find_neighbors();

  // Create an equation systems object, where ThinSystem
  // offers only the crucial functionality for solving a
  // system.  Use ThinSystem when you want the sleekest
  // system possible.
  EquationSystems equation_systems (mesh);

  // Declare the system and its variables.
  // Create a system named "Wave".  This can
  // be a simple, steady system
  equation_systems.add_system<LinearImplicitSystem> ("Wave");

  // Create an FEType describing the approximation
  // characteristics of the InfFE object.  Note that
  // the constructor automatically defaults to some
  // sensible values.  But use FIRST order
  // approximation.
  FEType fe_type(FIRST);

  // Add the variable "p" to "Wave".  Note that there exist
  // various approaches in adding variables.  In example 3,
  // add_variable took the order of approximation and used
  // default values for the FEFamily, while here the FEType
  // is used.
  equation_systems.get_system("Wave").add_variable("p", fe_type);

  // Give the system a pointer to the matrix assembly
  // function.
  equation_systems.get_system("Wave").attach_assemble_function (assemble_wave);

  // Set the speed of sound and fluid density
  // as EquationSystems parameter,
  // so that assemble_wave() can access it.
  equation_systems.parameters.set<Real>("speed")         = 1.;
  equation_systems.parameters.set<Real>("fluid density") = 1.;

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

#ifdef LIBMESH_ENABLE_AMR
  // Do uniform refinement if requested
  const unsigned int nr = input("nr", 0);
  if (nr)
    {
      MeshRefinement mesh_refinement(mesh);
      mesh_refinement.uniformly_refine(nr);
      equation_systems.reinit();
      equation_systems.print_info();
    }
#endif

  // Print and solve the refined sysem
  equation_systems.get_system("Wave").solve();

  libMesh::out << "Wave system solved" << std::endl;

  // Write the whole EquationSystems object to file.
  // For infinite elements, the concept of nodal_soln()
  // is not applicable. Therefore, writing the mesh in
  // some format @e always gives all-zero results at
  // the nodes of the infinite elements.  Instead,
  // use the FEInterface::compute_data() methods to
  // determine physically correct results within an
  // infinite element.
  equation_systems.write ("eqn_sys.dat", WRITE);

  libMesh::out << "eqn_sys.dat written" << std::endl;

  // All done.
  return 0;

#endif // else part of ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
}