Example 19 - Solving the 2D Young Laplace Problem using nonlinear solvers

This example shows how to use the NonlinearImplicitSystem class to efficiently solve nonlinear problems in parallel.

In nonlinear systems, we aim at finding x that satisfy R(x) = 0. In nonlinear finite element analysis, the residual is typically of the form R(x) = K(x)*x - f, with K(x) the system matrix and f the "right-hand-side". The NonlinearImplicitSystem class expects two callback functions to compute the residual R and its Jacobian for the Newton iterations. Here, we just approximate the true Jacobian by K(x).

You can turn on preconditining of the matrix free system using the jacobian by passing "-pre" on the command line. Currently this only work with Petsc so this isn't used by using "make run"

This example also runs with the experimental Trilinos NOX solvers by specifying the --use-trilinos command line argument.

C++ include files that we need

        #include <iostream>
        #include <algorithm>
        #include <cmath>

Various include files needed for the mesh & solver functionality.

        #include "libmesh.h"
        #include "mesh.h"
        #include "mesh_refinement.h"
        #include "gmv_io.h"
        #include "equation_systems.h"
        #include "fe.h"
        #include "quadrature_gauss.h"
        #include "dof_map.h"
        #include "sparse_matrix.h"
        #include "numeric_vector.h"
        #include "dense_matrix.h"
        #include "dense_vector.h"
        #include "elem.h"
        #include "string_to_enum.h"
        #include "getpot.h"

The nonlinear solver and system we will be using

        #include "nonlinear_solver.h"
        #include "nonlinear_implicit_system.h"

Necessary for programmatically setting petsc options

        #ifdef LIBMESH_HAVE_PETSC
        #include <petsc.h>
        #endif

A reference to our equation system

        EquationSystems *_equation_system = NULL;

Let-s define the physical parameters of the equation

        const Real kappa = 1.;
        const Real sigma = 0.2;

This function computes the Jacobian K(x)

        void compute_jacobian (const NumericVector<Number>& soln,
        		       SparseMatrix<Number>&  jacobian)
        {

Get a reference to the equation system.

          EquationSystems &es = *_equation_system;

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();

Get a reference to the NonlinearImplicitSystem we are solving

          NonlinearImplicitSystem& system = 
            es.get_system<NonlinearImplicitSystem>("Laplace-Young");

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. We will talk more about the \p 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. Since the \p FEBase::build() member dynamically creates memory we will store the object as an \p AutoPtr. 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);

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 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 Jacobian element matrix. Following basic finite element terminology we will denote these "Ke". More detail is in example 3.

          DenseMatrix<Number> Ke;

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 Jacobian contribution.

          MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
          const MeshBase::const_element_iterator end_el = mesh.active_local_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 Jacobian 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());

Now we will build the element Jacobian. This involves a double loop to integrate the test funcions (i) against the trial functions (j). Note that the Jacobian depends on the current solution x, which we access using the soln vector.

              for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        	{
        	  RealGradient grad_u = 0;
            
        	  for (unsigned int i=0; i<phi.size(); i++)
        	    grad_u += dphi[i][qp]*soln(dof_indices[i]);
        	  
        	  const Real K = 1./std::sqrt(1. + grad_u*grad_u);
        	  
        	  for (unsigned int i=0; i<phi.size(); i++)
        	    for (unsigned int j=0; j<phi.size(); j++)
        	      Ke(i,j) += JxW[qp]*(
        				  K*(dphi[i][qp]*dphi[j][qp]) +
        				  kappa*phi[i][qp]*phi[j][qp]
        				  );
        	}
              
              dof_map.constrain_element_matrix (Ke, dof_indices);

Add the element matrix to the system Jacobian.

              jacobian.add_matrix (Ke, dof_indices);
            }

That's it.

Here we compute the residual R(x) = K(x)*x - f. The current solution x is passed in the soln vector

        void compute_residual (const NumericVector<Number>& soln,
        		       NumericVector<Number>& residual)
        {
          EquationSystems &es = *_equation_system;

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();
          libmesh_assert (dim == 2);

Get a reference to the NonlinearImplicitSystem we are solving

          NonlinearImplicitSystem& system = 
            es.get_system<NonlinearImplicitSystem>("Laplace-Young");

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

          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 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 resdual contributions

          DenseVector<Number> Re;

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 residual.

          residual.zero();
        
          MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
          const MeshBase::const_element_iterator end_el = mesh.active_local_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);

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).

              Re.resize (dof_indices.size());

Now we will build the residual. This involves the construction of the matrix K and multiplication of it with the current solution x. We rearrange this into two loops: In the first, we calculate only the contribution of K_ij*x_j which is independent of the row i. In the second loops, we multiply with the row-dependent part and add it to the element residual.

              for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        	{
        	  Real u = 0;
        	  RealGradient grad_u = 0;
        	  
        	  for (unsigned int j=0; j<phi.size(); j++)
        	    {
        	      u      += phi[j][qp]*soln(dof_indices[j]);
        	      grad_u += dphi[j][qp]*soln(dof_indices[j]);
        	    }
        	  
        	  const Real K = 1./std::sqrt(1. + grad_u*grad_u);
        	  
        	  for (unsigned int i=0; i<phi.size(); i++)
        	    Re(i) += JxW[qp]*(
        			      K*(dphi[i][qp]*grad_u) +
        			      kappa*phi[i][qp]*u
        			      );
        	}

At this point the interior element integration has been completed. However, we have not yet addressed boundary conditions.

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

Loop over the face quadrature points for integration.

                    for (unsigned int qp=0; qp<qface.n_points(); qp++)
        	      {

This is the right-hand-side contribution (f), which has to be subtracted from the current residual

                        for (unsigned int i=0; i<phi_face.size(); i++)
        		  Re(i) -= JxW_face[qp]*sigma*phi_face[i][qp];
        	      } 
        	  }
              
              dof_map.constrain_element_vector (Re, dof_indices);
              residual.add_vector (Re, dof_indices);
            }

That's it.

Begin the main program.

        int main (int argc, char** argv)
        {

Initialize libMesh and any dependent libaries, like in example 2.

          LibMeshInit init (argc, argv);

Braces are used to force object scope, like in example 2

Create a GetPot object to parse the command line

            GetPot command_line (argc, argv);

Check for proper calling arguments.

            if (argc < 3)
              {
                if (libMesh::processor_id() == 0)
        	  std::cerr << "Usage:\n"
        		    <<"\t " << argv[0] << " -r 2"
        		    << std::endl;

This handy function will print the file name, line number, and then abort. Currrently the library does not use C++ exception handling.

                error();
              }

Brief message to the user regarding the program name and command line arguments.

            else 
              {
        	std::cout << "Running " << argv[0];
        	
        	for (int i=1; i<argc; i++)
        	  std::cout << " " << argv[i];
        	
        	std::cout << std::endl << std::endl;
              }

The dimension of our problem

            const unsigned int dim = 2;

Read number of refinements

            int nr = 2;
            if ( command_line.search(1, "-r") )
              nr = command_line.next(nr);

Read FE order from command line

            std::string order = "FIRST"; 
            if ( command_line.search(2, "-Order", "-o") )
              order = command_line.next(order);

Read FE Family from command line

            std::string family = "LAGRANGE"; 
            if ( command_line.search(2, "-FEFamily", "-f") )
              family = command_line.next(family);

Cannot use dicontinuous basis.

            if ((family == "MONOMIAL") || (family == "XYZ"))
              {
        	std::cout << "ex19 currently requires a C^0 (or higher) FE basis." << std::endl;
        	error();
              }
        
            if ( command_line.search(1, "-pre") )
              {
        #ifdef LIBMESH_HAVE_PETSC

Use the jacobian for preconditioning.

                PetscOptionsSetValue("-snes_mf_operator",PETSC_NULL);
        #else
                std::cerr<<"Must be using PetsC to use jacobian based preconditioning"<<std::endl;

returning zero so that "make run" won't fail if we ever enable this capability there.

                return 0;
        #endif //LIBMESH_HAVE_PETSC
              }

Create a mesh with user-defined dimension.

            Mesh mesh (dim);    
            mesh.read ("lshaped.xda");
        
            if (order != "FIRST")
              mesh.all_second_order();
        
            MeshRefinement(mesh).uniformly_refine(nr);

Print information about the mesh to the screen.

            mesh.print_info();

Create an equation systems object.

            EquationSystems equation_systems (mesh);
            _equation_system = &equation_systems;

Declare the system and its variables.

Creates a system named "Laplace-Young"

              NonlinearImplicitSystem& system =
        	equation_systems.add_system<NonlinearImplicitSystem> ("Laplace-Young");

Here we specify the tolerance for the nonlinear solver and the maximum of nonlinear iterations.

              equation_systems.parameters.set<Real>         ("nonlinear solver tolerance")          = 1.e-12;
              equation_systems.parameters.set<unsigned int> ("nonlinear solver maximum iterations") = 50;

Adds the variable "u" to "Laplace-Young". "u" will be approximated using second-order approximation.

              system.add_variable("u",
        			  Utility::string_to_enum<Order>   (order),
        			  Utility::string_to_enum<FEFamily>(family));

Give the system a pointer to the functions that update the residual and Jacobian.

              system.nonlinear_solver->residual = compute_residual;
              system.nonlinear_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();
            }

Solve the system "Laplace-Young"

            equation_systems.get_system("Laplace-Young").solve();

After solving the system write the solution

            GMVIO (mesh).write_equation_systems ("out.gmv", 
        					 equation_systems);
          }

All done.

          return 0; 
        }

The program without comments:

 
   
  
  #include <iostream>
  #include <algorithm>
  #include <cmath>
  
  #include "libmesh.h"
  #include "mesh.h"
  #include "mesh_refinement.h"
  #include "gmv_io.h"
  #include "equation_systems.h"
  #include "fe.h"
  #include "quadrature_gauss.h"
  #include "dof_map.h"
  #include "sparse_matrix.h"
  #include "numeric_vector.h"
  #include "dense_matrix.h"
  #include "dense_vector.h"
  #include "elem.h"
  #include "string_to_enum.h"
  #include "getpot.h"
  
  #include "nonlinear_solver.h"
  #include "nonlinear_implicit_system.h"
  
  #ifdef LIBMESH_HAVE_PETSC
  #include <petsc.h>
  #endif
  
  EquationSystems *_equation_system = NULL;
  
  const Real kappa = 1.;
  const Real sigma = 0.2;
  
  
  void compute_jacobian (const NumericVector<Number>& soln,
  		       SparseMatrix<Number>&  jacobian)
  {
    EquationSystems &es = *_equation_system;
  
    const MeshBase& mesh = es.get_mesh();
  
    const unsigned int dim = mesh.mesh_dimension();
  
    NonlinearImplicitSystem& system = 
      es.get_system<NonlinearImplicitSystem>("Laplace-Young");
    
    const DofMap& dof_map = system.get_dof_map();
  
    FEType fe_type = dof_map.variable_type(0);
  
    AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
    
    QGauss qrule (dim, FIFTH);
  
    fe->attach_quadrature_rule (&qrule);
  
    const std::vector<Real>& JxW = fe->get_JxW();
  
    const std::vector<std::vector<Real> >& phi = fe->get_phi();
    
    const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
  
    DenseMatrix<Number> Ke;
  
    std::vector<unsigned int> dof_indices;
  
    MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
    const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
  
    for ( ; el != end_el; ++el)
      {
        const Elem* elem = *el;
  
        dof_map.dof_indices (elem, dof_indices);
  
        fe->reinit (elem);
  
        Ke.resize (dof_indices.size(),
  		 dof_indices.size());
             
        for (unsigned int qp=0; qp<qrule.n_points(); qp++)
  	{
  	  RealGradient grad_u = 0;
      
  	  for (unsigned int i=0; i<phi.size(); i++)
  	    grad_u += dphi[i][qp]*soln(dof_indices[i]);
  	  
  	  const Real K = 1./std::sqrt(1. + grad_u*grad_u);
  	  
  	  for (unsigned int i=0; i<phi.size(); i++)
  	    for (unsigned int j=0; j<phi.size(); j++)
  	      Ke(i,j) += JxW[qp]*(
  				  K*(dphi[i][qp]*dphi[j][qp]) +
  				  kappa*phi[i][qp]*phi[j][qp]
  				  );
  	}
        
        dof_map.constrain_element_matrix (Ke, dof_indices);
        
        jacobian.add_matrix (Ke, dof_indices);
      }
  
  }
  
  
  void compute_residual (const NumericVector<Number>& soln,
  		       NumericVector<Number>& residual)
  {
    EquationSystems &es = *_equation_system;
  
    const MeshBase& mesh = es.get_mesh();
  
    const unsigned int dim = mesh.mesh_dimension();
    libmesh_assert (dim == 2);
  
    NonlinearImplicitSystem& system = 
      es.get_system<NonlinearImplicitSystem>("Laplace-Young");
    
    const DofMap& dof_map = system.get_dof_map();
  
    FEType fe_type = dof_map.variable_type(0);
  
    AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
    
    QGauss qrule (dim, FIFTH);
  
    fe->attach_quadrature_rule (&qrule);
  
    AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
  	      
    QGauss qface(dim-1, FIFTH);
    
    fe_face->attach_quadrature_rule (&qface);
  
    const std::vector<Real>& JxW = fe->get_JxW();
  
    const std::vector<std::vector<Real> >& phi = fe->get_phi();
    
    const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
  
    DenseVector<Number> Re;
  
    std::vector<unsigned int> dof_indices;
  
    residual.zero();
  
    MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
    const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
  
    for ( ; el != end_el; ++el)
      {
        const Elem* elem = *el;
  
        dof_map.dof_indices (elem, dof_indices);
  
        fe->reinit (elem);
  
        Re.resize (dof_indices.size());
        
  
        for (unsigned int qp=0; qp<qrule.n_points(); qp++)
  	{
  	  Real u = 0;
  	  RealGradient grad_u = 0;
  	  
  	  for (unsigned int j=0; j<phi.size(); j++)
  	    {
  	      u      += phi[j][qp]*soln(dof_indices[j]);
  	      grad_u += dphi[j][qp]*soln(dof_indices[j]);
  	    }
  	  
  	  const Real K = 1./std::sqrt(1. + grad_u*grad_u);
  	  
  	  for (unsigned int i=0; i<phi.size(); i++)
  	    Re(i) += JxW[qp]*(
  			      K*(dphi[i][qp]*grad_u) +
  			      kappa*phi[i][qp]*u
  			      );
  	}
  
        
        for (unsigned int side=0; side<elem->n_sides(); side++)
  	if (elem->neighbor(side) == NULL)
  	  {
  	    const std::vector<std::vector<Real> >&  phi_face = fe_face->get_phi();
  
  	    const std::vector<Real>& JxW_face = fe_face->get_JxW();
  
  	    fe_face->reinit(elem, side);
  
  	    for (unsigned int qp=0; qp<qface.n_points(); qp++)
  	      {
  		for (unsigned int i=0; i<phi_face.size(); i++)
  		  Re(i) -= JxW_face[qp]*sigma*phi_face[i][qp];
  	      } 
  	  }
        
        dof_map.constrain_element_vector (Re, dof_indices);
        residual.add_vector (Re, dof_indices);
      }
  
  }
  
  
  
  int main (int argc, char** argv)
  {
    LibMeshInit init (argc, argv);
   
    {
      GetPot command_line (argc, argv);
      
      if (argc < 3)
        {
          if (libMesh::processor_id() == 0)
  	  std::cerr << "Usage:\n"
  		    <<"\t " << argv[0] << " -r 2"
  		    << std::endl;
  
  	error();
        }
      
      else 
        {
  	std::cout << "Running " << argv[0];
  	
  	for (int i=1; i<argc; i++)
  	  std::cout << " " << argv[i];
  	
  	std::cout << std::endl << std::endl;
        }
      
  
      const unsigned int dim = 2;
      
      int nr = 2;
      if ( command_line.search(1, "-r") )
        nr = command_line.next(nr);
      
      std::string order = "FIRST"; 
      if ( command_line.search(2, "-Order", "-o") )
        order = command_line.next(order);
  
      std::string family = "LAGRANGE"; 
      if ( command_line.search(2, "-FEFamily", "-f") )
        family = command_line.next(family);
      
      if ((family == "MONOMIAL") || (family == "XYZ"))
        {
  	std::cout << "ex19 currently requires a C^0 (or higher) FE basis." << std::endl;
  	error();
        }
  
      if ( command_line.search(1, "-pre") )
        {
  #ifdef LIBMESH_HAVE_PETSC
          PetscOptionsSetValue("-snes_mf_operator",PETSC_NULL);
  #else
          std::cerr<<"Must be using PetsC to use jacobian based preconditioning"<<std::endl;
  
          return 0;
  #endif //LIBMESH_HAVE_PETSC
        }  
        
      Mesh mesh (dim);    
      mesh.read ("lshaped.xda");
  
      if (order != "FIRST")
        mesh.all_second_order();
  
      MeshRefinement(mesh).uniformly_refine(nr);
  
      mesh.print_info();    
      
      EquationSystems equation_systems (mesh);
      _equation_system = &equation_systems;
      
      {
        NonlinearImplicitSystem& system =
  	equation_systems.add_system<NonlinearImplicitSystem> ("Laplace-Young");
  
  
        equation_systems.parameters.set<Real>         ("nonlinear solver tolerance")          = 1.e-12;
        equation_systems.parameters.set<unsigned int> ("nonlinear solver maximum iterations") = 50;
  
        
        system.add_variable("u",
  			  Utility::string_to_enum<Order>   (order),
  			  Utility::string_to_enum<FEFamily>(family));
  
        
        system.nonlinear_solver->residual = compute_residual;
        system.nonlinear_solver->jacobian = compute_jacobian;
  
        equation_systems.init();
  
        equation_systems.print_info();
      }
      
      equation_systems.get_system("Laplace-Young").solve();
  
      GMVIO (mesh).write_equation_systems ("out.gmv", 
  					 equation_systems);
    }
    
    return 0; 
  }

The console output of the program:

***************************************************************
* Running Example  ./ex19-dbg
***************************************************************
 
Running ./ex19-dbg -r 3 -o FIRST

 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=225
    n_local_nodes()=225
  n_elem()=255
    n_local_elem()=255
    n_active_elem()=192
  n_subdomains()=1
  n_processors()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System "Laplace-Young"
    Type "NonlinearImplicit"
    Variables="u" 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="FIRST", "THIRD" 
    n_dofs()=225
    n_local_dofs()=225
    n_constrained_dofs()=0
    n_vectors()=1

  NL step 0, |residual|_2 = 2.000000e-01
  NL step 1, |residual|_2 = 4.432961e-03
  NL step 2, |residual|_2 = 2.163781e-04
  NL step 3, |residual|_2 = 1.157690e-05
  NL step 4, |residual|_2 = 6.567452e-07
  NL step 5, |residual|_2 = 3.849499e-08
  NL step 6, |residual|_2 = 2.293601e-09

 ---------------------------------------------------------------------------- 
| Reference count information                                                |
 ---------------------------------------------------------------------------- 
| 12SparseMatrixIdE reference count information:
|  Creations:    13
|  Destructions: 13
| 13NumericVectorIdE reference count information:
|  Creations:    23
|  Destructions: 23
| 15NonlinearSolverIdE reference count information:
|  Creations:    1
|  Destructions: 1
| 4Elem reference count information:
|  Creations:    1615
|  Destructions: 1615
| 4Node reference count information:
|  Creations:    225
|  Destructions: 225
| 5QBase reference count information:
|  Creations:    33
|  Destructions: 33
| 6DofMap reference count information:
|  Creations:    1
|  Destructions: 1
| 6FEBase reference count information:
|  Creations:    21
|  Destructions: 21
| 6System reference count information:
|  Creations:    1
|  Destructions: 1
| 9DofObject reference count information:
|  Creations:    1840
|  Destructions: 1840
| N10Parameters5ValueE reference count information:
|  Creations:    10
|  Destructions: 10
 ---------------------------------------------------------------------------- 

----------------------------------------------------------------------------------------------------------------
| Time:           Fri Oct 17 18:21:04 2008                                                                      |
| OS:             Darwin                                                                                        |
| HostName:       benjamin-kirks-macbook.local                                                                  |
| OS Release:     9.5.0                                                                                         |
| OS Version:     Darwin Kernel Version 9.5.0: Wed Sep  3 11:29:43 PDT 2008; root:xnu-1228.7.58~1/RELEASE_I386  |
| Machine:        i386                                                                                          |
| Username:       benkirk                                                                                       |
| Configuration:  ./configure run on Fri Oct 17 17:29:52 CDT 2008                                               |
----------------------------------------------------------------------------------------------------------------
 -------------------------------------------------------------------------------
| libMesh Performance: Alive time=0.633672, Active time=0.503336                |
 -------------------------------------------------------------------------------
| Event                          nCalls    Total       Avg         Percent of   |
|                                          Time        Time        Active Time  |
|-------------------------------------------------------------------------------|
|                                                                               |
|                                                                               |
| DofMap                                                                        |
|   add_neighbors_to_send_list() 1         0.0021      0.002114    0.42         |
|   compute_sparsity()           1         0.0643      0.064275    12.77        |
|   create_dof_constraints()     1         0.0018      0.001836    0.36         |
|   distribute_dofs()            1         0.0030      0.003020    0.60         |
|   dof_indices()                2880      0.0166      0.000006    3.29         |
|   reinit()                     1         0.0035      0.003548    0.70         |
|   sort_send_list()             1         0.0015      0.001458    0.29         |
|                                                                               |
| FE                                                                            |
|   compute_affine_map()         2944      0.0562      0.000019    11.17        |
|   compute_face_map()           448       0.0146      0.000033    2.90         |
|   compute_shape_functions()    2944      0.0653      0.000022    12.97        |
|   init_face_shape_functions()  161       0.0016      0.000010    0.31         |
|   init_shape_functions()       461       0.0159      0.000034    3.16         |
|   inverse_map()                2688      0.0399      0.000015    7.93         |
|                                                                               |
| GMVIO                                                                         |
|   write_nodal_data()           1         0.0032      0.003171    0.63         |
|                                                                               |
| LocationMap                                                                   |
|   find()                       756       0.0032      0.000004    0.64         |
|   init()                       3         0.0010      0.000332    0.20         |
|                                                                               |
| Mesh                                                                          |
|   find_neighbors()             2         0.0251      0.012527    4.98         |
|   renumber_nodes_and_elem()    2         0.0014      0.000716    0.28         |
|                                                                               |
| MeshRefinement                                                                |
|   _refine_elements()           3         0.0077      0.002561    1.53         |
|   add_point()                  756       0.0046      0.000006    0.91         |
|                                                                               |
| Parallel                                                                      |
|   allgather()                  1         0.0000      0.000008    0.00         |
|                                                                               |
| Partitioner                                                                   |
|   single_partition()           2         0.0016      0.000804    0.32         |
|                                                                               |
| System                                                                        |
|   assemble()                   1         0.0000      0.000003    0.00         |
|   solve()                      1         0.1693      0.169319    33.64        |
 -------------------------------------------------------------------------------
| Totals:                        14060     0.5033                  100.00       |
 -------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example  ./ex19-dbg
***************************************************************