// 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".  These datatypes are templated on  //  Number, which allows the same code to work for real  // or complex numbers.  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.  //  // Element iterators are a nice way to iterate through  // all the elements, or all the elements that have some property.  // There are many types of element iterators, but here we will  // use the most basic type, the  const_elem_iterator.  The iterator  //  el will iterate from the first to the last element.  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!//   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();  // Loop over the elements.  Note that  ++el is preferred to  // el++ since the latter requires an unnecessary temporary  // object.  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 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).      // The  DenseMatrix::resize() and the  DenseVector::resize()      // members will automatically zero out the matrix  and vector.      Ke.resize (dof_indices.size(),                 dof_indices.size());      Fe.resize (dof_indices.size());      // 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]);              }                    // This is the end of the matrix summation loop          // 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.          {            const Real x = q_point[qp](0);            const Real y = q_point[qp](1);            const Real eps = 1.e-3;                        // "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 is            //            // u_xx + u_yy = (u(i,j-1) + u(i,j+1) +            //                u(i-1,j) + u(i+1,j) +            //                -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 fxy = -(exact_solution(x,y-eps) +                               exact_solution(x,y+eps) +                               exact_solution(x-eps,y) +                               exact_solution(x+eps,y) -                               4.*exact_solution(x,y))/eps/eps;                        for (unsigned int i=0; i<phi.size(); i++)              Fe(i) += JxW[qp]*fxy*phi[i][qp];          }         }             // We have now reached the end of the RHS summation,      // and the end of quadrature point loop, so      // 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.      //      // There are several ways Dirichlet boundary conditions      // can be imposed.  A simple approach, which works for      // interpolary bases like the standard Lagrange polynomials,      // is to assign function values to the      // degrees of freedom living on the domain boundary. This      // works well for interpolary bases, but is more difficult      // when non-interpolary (e.g Legendre or Hierarchic) bases      // are used.      //      // Dirichlet boundary conditions can also be imposed with a      // "penalty" method.  In this case essentially the L2 projection      // of the boundary values are added to the matrix. The      // projection is multiplied by some large factor so that, in      // floating point arithmetic, the existing (smaller) entries      // in the matrix and right-hand-side are effectively ignored.      //      // This amounts to adding a term of the form (in latex notation)      //      // \frac{1}{\epsilon} \int_{\delta \Omega} \phi_i \phi_j = \frac{1}{\epsilon} \int_{\delta \Omega} u \phi_i      //      // where      //      // \frac{1}{\epsilon} is the penalty parameter, defined such that \epsilon << 1      {        // The following loop is 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 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);                  // The penalty value.  \frac{1}{\epsilon}                  // in the discussion above.                  const Real penalty = 1.e10;                  // The boundary value.                  const Real value = exact_solution(xf, yf);                                    // 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];                }             }      }            // We have now finished the quadrature point loop,      // and have therefore applied all the boundary conditions.      //      // 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  SparseMatrix::add_matrix()      // and  NumericVector::add_vector() members do this for us.      system.matrix->add_matrix (Ke, dof_indices);      system.rhs->add_vector    (Fe, dof_indices);    }    // All done!}