// 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.  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 5th order Gauss quadrature rule for numerical integration.  QGauss qrule (dim, FIFTH);  // Tell the finite element object to use our quadrature rule.  fe->attach_quadrature_rule (&qrule);  // 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, FIFTH);    // Tell the finte element object to use our  // quadrature rule.  fe_face->attach_quadrature_rule (&qface);  // Here we define some references to cell-specific data that  // will be used to assemble the linear system.  // We begin with the element Jacobian * quadrature weight at each  // integration point.     const std::vector<Real>& JxW = fe->get_JxW();  // The physical XY locations of the quadrature points on the element.  // These might be useful for evaluating spatially varying material  // properties at the quadrature points.  const std::vector<Point>& q_point = fe->get_xyz();  // 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();  // Define data structures to contain the element matrix  // and right-hand-side vector contribution.  Following  // basic finite element terminology we will denote these  // "Ke" and "Fe". More detail is in example 3.  DenseMatrix<Number> Ke;  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.  See example 3 for a discussion of the  // element iterators.  Here we use the \p const_local_elem_iterator  // to indicate we only want to loop over elements that are assigned  // to the local processor.  This allows each processor to compute  // its components of the global matrix.  //  // "PARALLEL CHANGE"  MeshBase::const_element_iterator       el     = mesh.local_elements_begin();  const MeshBase::const_element_iterator end_el = mesh.local_elements_end();  for ( ; el != end_el; ++el)    {      // Start logging the shape function initialization.      // This is done through a simple function call with      // the name of the event to log.      perf_log.push("elem init");            // 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 matrix and right-hand side 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 a      // triangle, now we are on a quadrilateral).      Ke.resize (dof_indices.size(),                 dof_indices.size());      Fe.resize (dof_indices.size());      // Stop logging the shape function initialization.      // If you forget to stop logging an event the PerfLog      // object will probably catch the error and abort.      perf_log.pop("elem init");            // Now we will build the element matrix.  This involves      // a double loop to integrate the test funcions (i) against      // the trial functions (j).      //      // We have split the numeric integration into two loops      // so that we can log the matrix and right-hand-side      // computation seperately.      //      // Now start logging the element matrix computation      perf_log.push ("Ke");      for (unsigned int qp=0; qp<qrule.n_points(); qp++)        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]);                  // Stop logging the matrix computation      perf_log.pop ("Ke");      // Now we build the element right-hand-side contribution.      // This involves a single loop in which we integrate the      // "forcing function" in the PDE against the test functions.      //      // Start logging the right-hand-side computation      perf_log.push ("Fe");            for (unsigned int qp=0; qp<qrule.n_points(); qp++)        {          // fxy is the forcing function for the Poisson equation.          // In this case we set fxy to be a finite difference          // Laplacian approximation to the (known) exact solution.          //          // We will use the second-order accurate FD Laplacian          // approximation, which in 2D on a structured grid is          //          // u_xx + u_yy = (u(i-1,j) + u(i+1,j) +          //                u(i,j-1) + u(i,j+1) +          //                -4*u(i,j))/h^2          //          // Since the value of the forcing function depends only          // on the location of the quadrature point (q_point[qp])          // we will compute it here, outside of the i-loop                    const Real x = q_point[qp](0);          const Real y = q_point[qp](1);          const Real z = q_point[qp](2);          const Real eps = 1.e-3;          const Real uxx = (exact_solution(x-eps,y,z) +                            exact_solution(x+eps,y,z) +                            -2.*exact_solution(x,y,z))/eps/eps;                        const Real uyy = (exact_solution(x,y-eps,z) +                            exact_solution(x,y+eps,z) +                            -2.*exact_solution(x,y,z))/eps/eps;                    const Real uzz = (exact_solution(x,y,z-eps) +                            exact_solution(x,y,z+eps) +                            -2.*exact_solution(x,y,z))/eps/eps;          Real fxy;          if(dim==1)          {            // In 1D, compute the rhs by differentiating the            // exact solution twice.            const Real pi = libMesh::pi;            fxy = (0.25*pi*pi)*sin(.5*pi*x);          }          else          {            fxy = - (uxx + uyy + ((dim==2) ? 0. : uzz));          }           // Add the RHS contribution          for (unsigned int i=0; i<phi.size(); i++)            Fe(i) += JxW[qp]*fxy*phi[i][qp];                  }            // Stop logging the right-hand-side computation      perf_log.pop ("Fe");      // At this point the interior element integration has      // been completed.  However, we have not yet addressed      // boundary conditions.  For this example we will only      // consider simple Dirichlet boundary conditions imposed      // via the penalty method. This is discussed at length in      // example 3.      {                // Start logging the boundary condition computation        perf_log.push ("BCs");        // The following loops over the sides of the element.        // If the element has no neighbor on a side then that        // side MUST live on a boundary of the domain.        for (unsigned int side=0; side<elem->n_sides(); side++)          if (elem->neighbor(side) == NULL)            {                          // The penalty value.  \frac{1}{\epsilon}              // in the discussion above.              const Real penalty = 1.e10;              // 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();              // The XYZ locations (in physical space) of the              // quadrature points on the face.  This is where              // we will interpolate the boundary value function.              const std::vector<Point >& qface_point = fe_face->get_xyz();              // Compute the shape function values on the element              // face.              fe_face->reinit(elem, side);              // Loop over the face quadrature points for integration.              for (unsigned int qp=0; qp<qface.n_points(); qp++)              {                // The location on the boundary of the current                // face quadrature point.                const Real xf = qface_point[qp](0);                const Real yf = qface_point[qp](1);                const Real zf = qface_point[qp](2);                // The boundary value.                const Real value = exact_solution(xf, yf, zf);                // Matrix contribution of the L2 projection.                 for (unsigned int i=0; i<phi_face.size(); i++)                  for (unsigned int j=0; j<phi_face.size(); j++)                    Ke(i,j) += JxW_face[qp]*penalty*phi_face[i][qp]*phi_face[j][qp];                // Right-hand-side contribution of the L2                // projection.                for (unsigned int i=0; i<phi_face.size(); i++)                  Fe(i) += JxW_face[qp]*penalty*value*phi_face[i][qp];              }             }                            // Stop logging the boundary condition computation        perf_log.pop ("BCs");      }             // The element matrix and right-hand-side are now built      // for this element.  Add them to the global matrix and      // right-hand-side vector.  The \p PetscMatrix::add_matrix()      // and \p PetscVector::add_vector() members do this for us.      // Start logging the insertion of the local (element)      // matrix and vector into the global matrix and vector      perf_log.push ("matrix insertion");            system.matrix->add_matrix (Ke, dof_indices);      system.rhs->add_vector    (Fe, dof_indices);      // Start logging the insertion of the local (element)      // matrix and vector into the global matrix and vector      perf_log.pop ("matrix insertion");    }  // That's it.  We don't need to do anything else to the  // PerfLog.  When it goes out of scope (at this function return)  // it will print its log to the screen. Pretty easy, huh?}