// Write nodal results to file.  The results can then      // be viewed with e.g. gnuplot (run gnuplot and type      // 'plot "pressure_node.res" with lines' in the command line)      res_out << t_time << "\t"              << global_displacement[dof_no]              << std::endl;    }#endif    // All done.    return 0;}// This function assembles the system matrix and right-hand-side// for our wave equation.void assemble_wave(EquationSystems& es,                   const std::string& system_name){    // It is a good idea to make sure we are assembling  // the proper system.  libmesh_assert (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 and fluid density  // to a local variable.  const Real speed = es.parameters.get<Real>("speed");  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  // \p FEBase::build() member dynamically creates memory we will  // store the object as an \p AutoPtr<FEBase>.  This can be thought  // of as a pointer that will clean up after itself.  AutoPtr<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 \p DofMap object for this system.  The \p 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<unsigned int> dof_indices;  // Now we will loop over all the elements in the mesh.  // We will compute the element matrix and right-hand-side  // contribution.//   const_elem_iterator           el (mesh.elements_begin());//   const const_elem_iterator end_el (mesh.elements_end());  MeshBase::const_element_iterator       el     = mesh.elements_begin();  const MeshBase::const_element_iterator end_el = mesh.elements_end();    for ( ; el != end_el; ++el)    {      // Store a pointer to the element we are currently      // working on.  This allows for nicer syntax later.      const Elem* elem = *el;      // 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 funcions (i) against          // the trial functions (j).          for (unsigned int i=0; i<phi.size(); i++)            for (unsigned int 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        for (unsigned int side=0; side<elem->n_sides(); side++)          if (!true)            // if (elem->neighbor(side) == NULL)            {              // Declare a special finite element object for              // boundary integration.              AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));                            // Boundary integration requires one quadraure rule,              // with dimensionality one less than the dimensionality              // of the element.              QGauss qface(dim-1, SECOND);                            // Tell the finte 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 (unsigned int 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(side) == NULL)                // In this example the Dirichlet boundary conditions will be         // imposed via panalty method after the        // system is assembled.              } // end boundary condition section                // 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    // All done!  return;}// This function applies the initial conditionsvoid apply_initial(EquationSystems& es,                   const std::string& system_name){  // Get a reference to our system, as before  NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);    // Numeric vectors for the pressure, velocity and acceleration  // values.  NumericVector<Number>&  pres_vec       = t_system.get_vector("displacement");  NumericVector<Number>&  vel_vec        = t_system.get_vector("velocity");  NumericVector<Number>&  acc_vec        = t_system.get_vector("acceleration");  // Assume our fluid to be at rest, which would  // also be the default conditions in class NewmarkSystem,  // but let us do it explicetly here.  pres_vec.zero();  vel_vec.zero();  acc_vec.zero();}// This function applies the Dirichlet boundary conditionsvoid fill_dirichlet_bc(EquationSystems& es,                       const std::string& system_name){  // It is a good idea to make sure we are assembling  // the proper system.  libmesh_assert (system_name == "Wave");  // Get a reference to our system, as before.  NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);  // Get writable references to the overall matrix and vector.  SparseMatrix<Number>&  matrix = *t_system.matrix;  NumericVector<Number>& rhs    = *t_system.rhs;  // Get a constant reference to the mesh object.  const MeshBase& mesh = es.get_mesh();  // Get \p libMesh's  \f$ \pi \f$  const Real pi = libMesh::pi;  // Ask the \p EquationSystems flag whether  // we should do this also for the matrix  const bool do_for_matrix =    es.parameters.get<bool>("Newmark set BC for Matrix");  // Ge the current time from \p EquationSystems  const Real current_time = es.parameters.get<Real>("time");  // Number of nodes in the mesh.  unsigned int n_nodes = mesh.n_nodes();  for (unsigned int n_cnt=0; n_cnt<n_nodes; n_cnt++)    {      // Get a reference to the current node.      const Node& curr_node = mesh.node(n_cnt);      // Check if Dirichlet BCs should be applied to this node.      // Use the \p TOLERANCE from \p mesh_common.h as tolerance.      // Here a pressure value is applied if the z-coord.      // is equal to 4, which corresponds to one end of the      // pipe-mesh in this directory.      const Real z_coo = 4.;      if (fabs(curr_node(2)-z_coo) < TOLERANCE)        {          // The global number of the respective degree of freedom.          unsigned int dn = curr_node.dof_number(0,0,0);          // The penalty parameter.          const Real penalty = 1.e10;          // Here we apply sinusoidal pressure values for 0<t<0.002          // at one end of the pipe-mesh.          Real p_value;          if (current_time < .002 )            p_value = sin(2*pi*current_time/.002);          else            p_value = .0;          // Now add the contributions to the matrix and the rhs.          rhs.add(dn, p_value*penalty);          // Add the panalty parameter to the global matrix          // if desired.          if (do_for_matrix)            matrix.add(dn, dn, penalty);        }    } // loop n_cnt}