The source file exact_solution.C with comments:

        #include <math.h>
        
Mesh library includes
        #include "libmesh/libmesh_common.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
        
        
        
        
        /**
         * This is the exact solution that
         * we are trying to obtain.  We will solve
         *
         * - (u_xx + u_yy) = f
         *
         * and take a finite difference approximation using this
         * function to get f.  This is the well-known "method of
         * manufactured solutions".
         */
        Real exact_solution (const Real x,
        		     const Real y,
        		     const Real z = 0.)
        {
          static const Real pi = acos(-1.);
        
          return cos(.5*pi*x)*sin(.5*pi*y)*cos(.5*pi*z);
        }



The source file introduction_ex3.C with comments:

Introduction Example 3 - Solving a Poisson Problem



This is the third example program. It builds on the second example program by showing how to solve a simple Poisson system. This example also introduces the notion of customized matrix assembly functions, working with an exact solution, and using element iterators. We will not comment on things that were already explained in the second example.

C++ include files that we need
        #include <iostream>
        #include <algorithm>
        #include <math.h>
        
Basic include files needed for the mesh functionality.
        #include "libmesh/libmesh.h"
        #include "libmesh/mesh.h"
        #include "libmesh/mesh_generation.h"
        #include "libmesh/vtk_io.h"
        #include "libmesh/linear_implicit_system.h"
        #include "libmesh/equation_systems.h"
        
Define the Finite Element object.
        #include "libmesh/fe.h"
        
Define Gauss quadrature rules.
        #include "libmesh/quadrature_gauss.h"
        
Define useful datatypes for finite element matrix and vector components.
        #include "libmesh/sparse_matrix.h"
        #include "libmesh/numeric_vector.h"
        #include "libmesh/dense_matrix.h"
        #include "libmesh/dense_vector.h"
        #include "libmesh/elem.h"
        
Define the DofMap, which handles degree of freedom indexing.
        #include "libmesh/dof_map.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
Function prototype. This is the function that will assemble the linear system for our Poisson problem. Note that the function will take the EquationSystems object and the name of the system we are assembling as input. From the EquationSystems object we have access to the Mesh and other objects we might need.
        void assemble_poisson(EquationSystems& es,
                              const std::string& system_name);
        
Function prototype for the exact solution.
        Real exact_solution (const Real x,
                             const Real y,
                             const Real z = 0.);
        
        int main (int argc, char** argv)
        {
Initialize libraries, like in example 2.
          LibMeshInit init (argc, argv);
        
Brief message to the user regarding the program name and command line arguments.
          std::cout << "Running " << argv[0];
        
          for (int i=1; i<argc; i++)
            std::cout << " " << argv[i];
        
          std::cout << std::endl << std::endl;
        
Skip this 2D example if libMesh was compiled as 1D-only.
          libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
        
Create a mesh, with dimension to be overridden later, distributed across the default MPI communicator.
          Mesh mesh(init.comm());
        
Use the MeshTools::Generation mesh generator to create a uniform 2D grid on the square [-1,1]^2. We instruct the mesh generator to build a mesh of 15x15 QUAD9 elements. Building QUAD9 elements instead of the default QUAD4's we used in example 2 allow us to use higher-order approximation.
          MeshTools::Generation::build_square (mesh,
                                               15, 15,
                                               -1., 1.,
                                               -1., 1.,
                                               QUAD9);
        
Print information about the mesh to the screen. Note that 5x5 QUAD9 elements actually has 11x11 nodes, so this mesh is significantly larger than the one in example 2.
          mesh.print_info();
        
Create an equation systems object.
          EquationSystems equation_systems (mesh);
        
Declare the Poisson system and its variables. The Poisson system is another example of a steady system.
          equation_systems.add_system<LinearImplicitSystem> ("Poisson");
        
Adds the variable "u" to "Poisson". "u" will be approximated using second-order approximation.
          equation_systems.get_system("Poisson").add_variable("u", SECOND);
        
Give the system a pointer to the matrix assembly function. This will be called when needed by the library.
          equation_systems.get_system("Poisson").attach_assemble_function (assemble_poisson);
        
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 "Poisson". Note that calling this member will assemble the linear system and invoke the default numerical solver. With PETSc the solver can be controlled from the command line. For example, you can invoke conjugate gradient with:

./introduction_ex3 -ksp_type cg

You can also get a nice X-window that monitors the solver convergence with:

./introduction-ex3 -ksp_xmonitor

if you linked against the appropriate X libraries when you built PETSc.
          equation_systems.get_system("Poisson").solve();
        
        #if defined(LIBMESH_HAVE_VTK) && !defined(LIBMESH_ENABLE_PARMESH)
        
After solving the system write the solution to a VTK-formatted plot file.
          VTKIO (mesh).write_equation_systems ("out.pvtu", equation_systems);
        
        #endif // #ifdef LIBMESH_HAVE_VTK
        
All done.
          return 0;
        }
        
        
        
We now define the matrix assembly function for the Poisson system. We need to first compute element matrices and right-hand sides, and then take into account the boundary conditions, which will be handled via a penalty method.
        void assemble_poisson(EquationSystems& es,
                              const std::string& system_name)
        {
        
It is a good idea to make sure we are assembling the proper system.
          libmesh_assert_equal_to (system_name, "Poisson");
        
        
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 LinearImplicitSystem we are solving
          LinearImplicitSystem& system = es.get_system<LinearImplicitSystem> ("Poisson");
        
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. We will talk more about the 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 FEBase::build() member dynamically creates memory we will store the object as an AutoPtr. This can be thought of as a pointer that will clean up after itself. Introduction Example 4 describes some advantages of AutoPtr's in the context of quadrature rules.
          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 finite 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.

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

Element iterators are a nice way to iterate through all the elements, or all the elements that have some property. The iterator el will iterate from the first to the last element on the local processor. 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! In case users later modify this program to include refinement, we will be safe and will only consider the active elements; hence we use a variant of the \p active_elem_iterator.
          MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
          const MeshBase::const_element_iterator end_el = mesh.active_local_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.

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_and_vector (Ke, Fe, dof_indices);
        
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!
        }



The source file exact_solution.C without comments:

 
  #include <math.h>
  
  #include "libmesh/libmesh_common.h"
  
  using namespace libMesh;
  
  
  
  
  
  /**
   * This is the exact solution that
   * we are trying to obtain.  We will solve
   *
   * - (u_xx + u_yy) = f
   *
   * and take a finite difference approximation using this
   * function to get f.  This is the well-known "method of
   * manufactured solutions".
   */
  Real exact_solution (const Real x,
  		     const Real y,
  		     const Real z = 0.)
  {
    static const Real pi = acos(-1.);
  
    return cos(.5*pi*x)*sin(.5*pi*y)*cos(.5*pi*z);
  }



The source file introduction_ex3.C without comments:

 
  
  #include <iostream>
  #include <algorithm>
  #include <math.h>
  
  #include "libmesh/libmesh.h"
  #include "libmesh/mesh.h"
  #include "libmesh/mesh_generation.h"
  #include "libmesh/vtk_io.h"
  #include "libmesh/linear_implicit_system.h"
  #include "libmesh/equation_systems.h"
  
  #include "libmesh/fe.h"
  
  #include "libmesh/quadrature_gauss.h"
  
  #include "libmesh/sparse_matrix.h"
  #include "libmesh/numeric_vector.h"
  #include "libmesh/dense_matrix.h"
  #include "libmesh/dense_vector.h"
  #include "libmesh/elem.h"
  
  #include "libmesh/dof_map.h"
  
  using namespace libMesh;
  
  void assemble_poisson(EquationSystems& es,
                        const std::string& system_name);
  
  Real exact_solution (const Real x,
                       const Real y,
                       const Real z = 0.);
  
  int main (int argc, char** argv)
  {
    LibMeshInit init (argc, argv);
  
    std::cout << "Running " << argv[0];
  
    for (int i=1; i<argc; i++)
      std::cout << " " << argv[i];
  
    std::cout << std::endl << std::endl;
  
    libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
  
    Mesh mesh(init.comm());
  
    MeshTools::Generation::build_square (mesh,
                                         15, 15,
                                         -1., 1.,
                                         -1., 1.,
                                         QUAD9);
  
    mesh.print_info();
  
    EquationSystems equation_systems (mesh);
  
    equation_systems.add_system<LinearImplicitSystem> ("Poisson");
  
    equation_systems.get_system("Poisson").add_variable("u", SECOND);
  
    equation_systems.get_system("Poisson").attach_assemble_function (assemble_poisson);
  
    equation_systems.init();
  
    equation_systems.print_info();
  
    equation_systems.get_system("Poisson").solve();
  
  #if defined(LIBMESH_HAVE_VTK) && !defined(LIBMESH_ENABLE_PARMESH)
  
    VTKIO (mesh).write_equation_systems ("out.pvtu", equation_systems);
  
  #endif // #ifdef LIBMESH_HAVE_VTK
  
    return 0;
  }
  
  
  
  void assemble_poisson(EquationSystems& es,
                        const std::string& system_name)
  {
  
    libmesh_assert_equal_to (system_name, "Poisson");
  
  
    const MeshBase& mesh = es.get_mesh();
  
    const unsigned int dim = mesh.mesh_dimension();
  
    LinearImplicitSystem& system = es.get_system<LinearImplicitSystem> ("Poisson");
  
    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<Point>& q_point = fe->get_xyz();
  
    const std::vector<std::vector<Real> >& phi = fe->get_phi();
  
    const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
  
    DenseMatrix<Number> Ke;
    DenseVector<Number> Fe;
  
  
    std::vector<dof_id_type> 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());
  
        Fe.resize (dof_indices.size());
  
        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]);
                }
  
            {
              const Real x = q_point[qp](0);
              const Real y = q_point[qp](1);
              const Real eps = 1.e-3;
  
  
              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];
            }
          }
  
        {
  
          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();
  
                const std::vector<Point >& qface_point = fe_face->get_xyz();
  
                fe_face->reinit(elem, side);
  
                for (unsigned int qp=0; qp<qface.n_points(); qp++)
                  {
  
                    const Real xf = qface_point[qp](0);
                    const Real yf = qface_point[qp](1);
  
                    const Real penalty = 1.e10;
  
                    const Real value = exact_solution(xf, yf);
  
                    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];
  
                    for (unsigned int i=0; i<phi_face.size(); i++)
                      Fe(i) += JxW_face[qp]*penalty*value*phi_face[i][qp];
                  }
              }
        }
  
  
        dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
  
        system.matrix->add_matrix (Ke, dof_indices);
        system.rhs->add_vector    (Fe, dof_indices);
      }
  
  }



The console output of the program:

make[4]: Entering directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/introduction/introduction_ex3'
***************************************************************
* Running Example introduction_ex3:
*  mpirun -np 4 example-devel  -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
 
Running /net/spark/workspace/roystgnr/libmesh/git/devel/examples/introduction/introduction_ex3/.libs/lt-example-devel -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc

 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=961
    n_local_nodes()=255
  n_elem()=225
    n_local_elem()=56
    n_active_elem()=225
  n_subdomains()=1
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System #0, "Poisson"
    Type "LinearImplicit"
    Variables="u" 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="SECOND", "THIRD" 
    n_dofs()=961
    n_local_dofs()=255
    n_constrained_dofs()=0
    n_local_constrained_dofs()=0
    n_vectors()=1
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 14.4849
      Average Off-Processor Bandwidth <= 1.01145
      Maximum  On-Processor Bandwidth <= 25
      Maximum Off-Processor Bandwidth <= 16
    DofMap Constraints
      Number of DoF Constraints = 0
      Number of Node Constraints = 0

*** Warning, This code is untested, experimental, or likely to see future API changes: ../src/mesh/vtk_io.C, line 358, compiled Apr 19 2013 at 11:42:41 ***

 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:44:33 2013                                                                          |
| OS:             Linux                                                                                             |
| HostName:       spark.ices.utexas.edu                                                                             |
| OS Release:     2.6.32-279.22.1.el6.x86_64                                                                        |
| OS Version:     #1 SMP Tue Feb 5 14:33:39 CST 2013                                                                |
| Machine:        x86_64                                                                                            |
| Username:       roystgnr                                                                                          |
| Configuration:  ../configure  '--enable-everything'                                                               |
|  'METHODS=devel'                                                                                                  |
|  '--prefix=/h2/roystgnr/libmesh-test'                                                                             |
|  'CXX=distcc /usr/bin/g++'                                                                                        |
|  'CC=distcc /usr/bin/gcc'                                                                                         |
|  'FC=distcc /usr/bin/gfortran'                                                                                    |
|  'F77=distcc /usr/bin/gfortran'                                                                                   |
|  'PETSC_DIR=/opt/apps/ossw/libraries/petsc/petsc-3.3-p2'                                                          |
|  'PETSC_ARCH=gcc-system-mkl-gf-10.3.12.361-mpich2-1.4.1p1-cxx-opt'                                                |
|  'SLEPC_DIR=/opt/apps/ossw/libraries/slepc/slepc-3.3-p2-petsc-3.3-p2-cxx-opt'                                     |
|  'TRILINOS_DIR=/opt/apps/ossw/libraries/trilinos/trilinos-10.12.2/sl6/gcc-system/mpich2-1.4.1p1/mkl-gf-10.3.12.361'|
|  'VTK_DIR=/opt/apps/ossw/libraries/vtk/vtk-5.10.0/sl6/gcc-system'                                                 |
|  'HDF5_DIR=/opt/apps/ossw/libraries/hdf5/hdf5-1.8.9/sl6/gcc-system'                                               |
 -------------------------------------------------------------------------------------------------------------------
 ----------------------------------------------------------------------------------------------------------------
| libMesh Performance: Alive time=0.051717, Active time=0.043031                                                 |
 ----------------------------------------------------------------------------------------------------------------
| Event                              nCalls    Total Time  Avg Time    Total Time  Avg Time    % of Active Time  |
|                                              w/o Sub     w/o Sub     With Sub    With Sub    w/o S    With S   |
|----------------------------------------------------------------------------------------------------------------|
|                                                                                                                |
|                                                                                                                |
| DofMap                                                                                                         |
|   add_neighbors_to_send_list()     1         0.0003      0.000335    0.0005      0.000510    0.78     1.19     |
|   build_sparsity()                 1         0.0003      0.000313    0.0013      0.001284    0.73     2.98     |
|   create_dof_constraints()         1         0.0001      0.000085    0.0001      0.000085    0.20     0.20     |
|   distribute_dofs()                1         0.0009      0.000911    0.0047      0.004691    2.12     10.90    |
|   dof_indices()                    199       0.0011      0.000006    0.0011      0.000006    2.64     2.64     |
|   prepare_send_list()              1         0.0000      0.000009    0.0000      0.000009    0.02     0.02     |
|   reinit()                         1         0.0012      0.001214    0.0012      0.001214    2.82     2.82     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          1         0.0002      0.000198    0.0012      0.001190    0.46     2.77     |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        71        0.0004      0.000006    0.0004      0.000006    0.91     0.91     |
|   init_shape_functions()           16        0.0000      0.000003    0.0000      0.000003    0.10     0.10     |
|   inverse_map()                    45        0.0002      0.000004    0.0002      0.000004    0.43     0.43     |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             71        0.0002      0.000003    0.0002      0.000003    0.56     0.56     |
|   compute_face_map()               15        0.0001      0.000008    0.0003      0.000020    0.27     0.70     |
|   init_face_shape_functions()      1         0.0000      0.000005    0.0000      0.000005    0.01     0.01     |
|   init_reference_to_physical_map() 16        0.0001      0.000009    0.0001      0.000009    0.33     0.33     |
|                                                                                                                |
| Mesh                                                                                                           |
|   find_neighbors()                 1         0.0004      0.000421    0.0008      0.000771    0.98     1.79     |
|   renumber_nodes_and_elem()        2         0.0001      0.000043    0.0001      0.000043    0.20     0.20     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   compute_hilbert_indices()        2         0.0012      0.000577    0.0012      0.000577    2.68     2.68     |
|   find_global_indices()            2         0.0002      0.000109    0.0029      0.001458    0.51     6.78     |
|   parallel_sort()                  2         0.0003      0.000155    0.0012      0.000615    0.72     2.86     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         1         0.0031      0.003082    0.0043      0.004342    7.16     10.09    |
|                                                                                                                |
| MeshTools::Generation                                                                                          |
|   build_cube()                     1         0.0003      0.000333    0.0003      0.000333    0.77     0.77     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      1         0.0017      0.001717    0.0033      0.003261    3.99     7.58     |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      9         0.0018      0.000202    0.0019      0.000207    4.23     4.33     |
|   max(bool)                        1         0.0000      0.000004    0.0000      0.000004    0.01     0.01     |
|   max(scalar)                      105       0.0005      0.000004    0.0005      0.000004    1.09     1.09     |
|   max(vector)                      24        0.0001      0.000005    0.0004      0.000017    0.31     0.96     |
|   min(bool)                        121       0.0004      0.000004    0.0004      0.000004    1.03     1.03     |
|   min(scalar)                      99        0.0058      0.000059    0.0058      0.000059    13.58    13.58    |
|   min(vector)                      24        0.0002      0.000008    0.0007      0.000030    0.43     1.69     |
|   probe()                          36        0.0006      0.000018    0.0006      0.000018    1.49     1.49     |
|   receive()                        36        0.0001      0.000002    0.0007      0.000020    0.21     1.70     |
|   send()                           36        0.0000      0.000001    0.0000      0.000001    0.12     0.12     |
|   send_receive()                   40        0.0002      0.000004    0.0010      0.000025    0.37     2.30     |
|   sum()                            20        0.0010      0.000051    0.0015      0.000076    2.39     3.53     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           36        0.0000      0.000001    0.0000      0.000001    0.09     0.09     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         1         0.0003      0.000265    0.0015      0.001455    0.62     3.38     |
|   set_parent_processor_ids()       1         0.0000      0.000044    0.0000      0.000044    0.10     0.10     |
|                                                                                                                |
| PetscLinearSolver                                                                                              |
|   solve()                          1         0.0176      0.017635    0.0176      0.017635    40.98    40.98    |
|                                                                                                                |
| System                                                                                                         |
|   assemble()                       1         0.0015      0.001531    0.0030      0.003023    3.56     7.03     |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            1044      0.0430                                          100.00            |
 ----------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example introduction_ex3:
*  mpirun -np 4 example-devel  -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
make[4]: Leaving directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/introduction/introduction_ex3'

Site Created By: libMesh Developers
Last modified: April 23 2013 04:19:30 UTC

Hosted By:
SourceForge.net Logo