The source file systems_of_equations_ex1.C with comments:

Systems Example 1 - Stokes Equations



This example shows how a simple, linear system of equations can be solved in parallel. The system of equations are the familiar Stokes equations for low-speed incompressible fluid flow.

C++ include files that we need
        #include <iostream>
        #include <algorithm>
        #include <math.h>
        
Basic include file needed for the mesh functionality.
        #include "libmesh/libmesh.h"
        #include "libmesh/mesh.h"
        #include "libmesh/mesh_generation.h"
        #include "libmesh/exodusII_io.h"
        #include "libmesh/equation_systems.h"
        #include "libmesh/fe.h"
        #include "libmesh/quadrature_gauss.h"
        #include "libmesh/dof_map.h"
        #include "libmesh/sparse_matrix.h"
        #include "libmesh/numeric_vector.h"
        #include "libmesh/dense_matrix.h"
        #include "libmesh/dense_vector.h"
        #include "libmesh/linear_implicit_system.h"
        
For systems of equations the \p DenseSubMatrix and \p DenseSubVector provide convenient ways for assembling the element matrix and vector on a component-by-component basis.
        #include "libmesh/dense_submatrix.h"
        #include "libmesh/dense_subvector.h"
        
The definition of a geometric element
        #include "libmesh/elem.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
Function prototype. This function will assemble the system matrix and right-hand-side.
        void assemble_stokes (EquationSystems& es,
                              const std::string& system_name);
        
The main program.
        int main (int argc, char** argv)
        {
Initialize libMesh.
          LibMeshInit init (argc, argv);
        
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 8x8 \p Quad9 elements. Building these higher-order elements allows us to use higher-order approximation, as in example 3.
          MeshTools::Generation::build_square (mesh,
                                               15, 15,
                                               0., 1.,
                                               0., 1.,
                                               QUAD9);
        
Print information about the mesh to the screen.
          mesh.print_info();
        
Create an equation systems object.
          EquationSystems equation_systems (mesh);
        
Declare the system and its variables. Create a transient system named "Convection-Diffusion"
          LinearImplicitSystem & system =
            equation_systems.add_system<LinearImplicitSystem> ("Stokes");
        
Add the variables "u" & "v" to "Stokes". They will be approximated using second-order approximation.
          system.add_variable ("u", SECOND);
          system.add_variable ("v", SECOND);
        
Add the variable "p" to "Stokes". This will be approximated with a first-order basis, providing an LBB-stable pressure-velocity pair.
          system.add_variable ("p", FIRST);
        
Give the system a pointer to the matrix assembly function.
          system.attach_assemble_function (assemble_stokes);
        
Initialize the data structures for the equation system.
          equation_systems.init ();
        
          equation_systems.parameters.set<unsigned int>("linear solver maximum iterations") = 250;
          equation_systems.parameters.set<Real>        ("linear solver tolerance") = TOLERANCE;
        
Prints information about the system to the screen.
          equation_systems.print_info();
        
Assemble & solve the linear system, then write the solution.
          equation_systems.get_system("Stokes").solve();
        
        #ifdef LIBMESH_HAVE_EXODUS_API
          ExodusII_IO(mesh).write_equation_systems ("out.e",
                                              equation_systems);
        #endif // #ifdef LIBMESH_HAVE_EXODUS_API
        
All done.
          return 0;
        }
        
        void assemble_stokes (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, "Stokes");
        
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 Convection-Diffusion system object.
          LinearImplicitSystem & system =
            es.get_system<LinearImplicitSystem> ("Stokes");
        
Numeric ids corresponding to each variable in the system
          const unsigned int u_var = system.variable_number ("u");
          const unsigned int v_var = system.variable_number ("v");
          const unsigned int p_var = system.variable_number ("p");
        
Get the Finite Element type for "u". Note this will be the same as the type for "v".
          FEType fe_vel_type = system.variable_type(u_var);
        
Get the Finite Element type for "p".
          FEType fe_pres_type = system.variable_type(p_var);
        
Build a Finite Element object of the specified type for the velocity variables.
          AutoPtr<FEBase> fe_vel  (FEBase::build(dim, fe_vel_type));
        
Build a Finite Element object of the specified type for the pressure variables.
          AutoPtr<FEBase> fe_pres (FEBase::build(dim, fe_pres_type));
        
A Gauss quadrature rule for numerical integration. Let the \p FEType object decide what order rule is appropriate.
          QGauss qrule (dim, fe_vel_type.default_quadrature_order());
        
Tell the finite element objects to use our quadrature rule.
          fe_vel->attach_quadrature_rule (&qrule);
          fe_pres->attach_quadrature_rule (&qrule);
        
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_vel->get_JxW();
        
The element shape function gradients for the velocity variables evaluated at the quadrature points.
          const std::vector<std::vector<RealGradient> >& dphi = fe_vel->get_dphi();
        
The element shape functions for the pressure variable evaluated at the quadrature points.
          const std::vector<std::vector<Real> >& psi = fe_pres->get_phi();
        
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();
        
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".
          DenseMatrix<Number> Ke;
          DenseVector<Number> Fe;
        
          DenseSubMatrix<Number>
            Kuu(Ke), Kuv(Ke), Kup(Ke),
            Kvu(Ke), Kvv(Ke), Kvp(Ke),
            Kpu(Ke), Kpv(Ke), Kpp(Ke);
        
          DenseSubVector<Number>
            Fu(Fe),
            Fv(Fe),
            Fp(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;
          std::vector<dof_id_type> dof_indices_u;
          std::vector<dof_id_type> dof_indices_v;
          std::vector<dof_id_type> dof_indices_p;
        
Now we will loop over all the elements in the mesh that live on the local processor. We will compute the element matrix and right-hand-side contribution. 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();
        
          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);
              dof_map.dof_indices (elem, dof_indices_u, u_var);
              dof_map.dof_indices (elem, dof_indices_v, v_var);
              dof_map.dof_indices (elem, dof_indices_p, p_var);
        
              const unsigned int n_dofs   = dof_indices.size();
              const unsigned int n_u_dofs = dof_indices_u.size();
              const unsigned int n_v_dofs = dof_indices_v.size();
              const unsigned int n_p_dofs = dof_indices_p.size();
        
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_vel->reinit  (elem);
              fe_pres->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 (n_dofs, n_dofs);
              Fe.resize (n_dofs);
        
Reposition the submatrices... The idea is this:

- - - - | Kuu Kuv Kup | | Fu | Ke = | Kvu Kvv Kvp |; Fe = | Fv | | Kpu Kpv Kpp | | Fp | - - - -

The \p DenseSubMatrix.repostition () member takes the (row_offset, column_offset, row_size, column_size).

Similarly, the \p DenseSubVector.reposition () member takes the (row_offset, row_size)
              Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
              Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
              Kup.reposition (u_var*n_u_dofs, p_var*n_u_dofs, n_u_dofs, n_p_dofs);
        
              Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
              Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
              Kvp.reposition (v_var*n_v_dofs, p_var*n_v_dofs, n_v_dofs, n_p_dofs);
        
              Kpu.reposition (p_var*n_u_dofs, u_var*n_u_dofs, n_p_dofs, n_u_dofs);
              Kpv.reposition (p_var*n_u_dofs, v_var*n_u_dofs, n_p_dofs, n_v_dofs);
              Kpp.reposition (p_var*n_u_dofs, p_var*n_u_dofs, n_p_dofs, n_p_dofs);
        
              Fu.reposition (u_var*n_u_dofs, n_u_dofs);
              Fv.reposition (v_var*n_u_dofs, n_v_dofs);
              Fp.reposition (p_var*n_u_dofs, n_p_dofs);
        
Now we will build the element matrix.
              for (unsigned int qp=0; qp<qrule.n_points(); qp++)
                {
Assemble the u-velocity row uu coupling
                  for (unsigned int i=0; i<n_u_dofs; i++)
                    for (unsigned int j=0; j<n_u_dofs; j++)
                      Kuu(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
        
up coupling
                  for (unsigned int i=0; i<n_u_dofs; i++)
                    for (unsigned int j=0; j<n_p_dofs; j++)
                      Kup(i,j) += -JxW[qp]*psi[j][qp]*dphi[i][qp](0);
        
        
Assemble the v-velocity row vv coupling
                  for (unsigned int i=0; i<n_v_dofs; i++)
                    for (unsigned int j=0; j<n_v_dofs; j++)
                      Kvv(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
        
vp coupling
                  for (unsigned int i=0; i<n_v_dofs; i++)
                    for (unsigned int j=0; j<n_p_dofs; j++)
                      Kvp(i,j) += -JxW[qp]*psi[j][qp]*dphi[i][qp](1);
        
        
Assemble the pressure row pu coupling
                  for (unsigned int i=0; i<n_p_dofs; i++)
                    for (unsigned int j=0; j<n_u_dofs; j++)
                      Kpu(i,j) += -JxW[qp]*psi[i][qp]*dphi[j][qp](0);
        
pv coupling
                  for (unsigned int i=0; i<n_p_dofs; i++)
                    for (unsigned int j=0; j<n_v_dofs; j++)
                      Kpv(i,j) += -JxW[qp]*psi[i][qp]*dphi[j][qp](1);
        
                } // end of the quadrature point qp-loop
        
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. The penalty method used here is equivalent (for Lagrange basis functions) to lumping the matrix resulting from the L2 projection penalty approach introduced in example 3.
              {
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 s=0; s<elem->n_sides(); s++)
                  if (elem->neighbor(s) == NULL)
                    {
                      AutoPtr<Elem> side (elem->build_side(s));
        
Loop over the nodes on the side.
                      for (unsigned int ns=0; ns<side->n_nodes(); ns++)
                        {
The location on the boundary of the current node.

const Real xf = side->point(ns)(0);
                          const Real yf = side->point(ns)(1);
        
The penalty value. \f$ \frac{1}{\epsilon \f$
                          const Real penalty = 1.e10;
        
The boundary values.

Set u = 1 on the top boundary, 0 everywhere else
                          const Real u_value = (yf > .99) ? 1. : 0.;
        
Set v = 0 everywhere
                          const Real v_value = 0.;
        
Find the node on the element matching this node on the side. That defined where in the element matrix the boundary condition will be applied.
                          for (unsigned int n=0; n<elem->n_nodes(); n++)
                            if (elem->node(n) == side->node(ns))
                              {
Matrix contribution.
                                Kuu(n,n) += penalty;
                                Kvv(n,n) += penalty;
        
Right-hand-side contribution.
                                Fu(n) += penalty*u_value;
                                Fv(n) += penalty*v_value;
                              }
                        } // end face node loop
                    } // end if (elem->neighbor(side) == NULL)
              } // end boundary condition section
        
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 \p NumericMatrix::add_matrix() and \p NumericVector::add_vector() members do this for us.
              system.matrix->add_matrix (Ke, dof_indices);
              system.rhs->add_vector    (Fe, dof_indices);
            } // end of element loop
        
That's it.
          return;
        }



The source file systems_of_equations_ex1.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/exodusII_io.h"
  #include "libmesh/equation_systems.h"
  #include "libmesh/fe.h"
  #include "libmesh/quadrature_gauss.h"
  #include "libmesh/dof_map.h"
  #include "libmesh/sparse_matrix.h"
  #include "libmesh/numeric_vector.h"
  #include "libmesh/dense_matrix.h"
  #include "libmesh/dense_vector.h"
  #include "libmesh/linear_implicit_system.h"
  
  #include "libmesh/dense_submatrix.h"
  #include "libmesh/dense_subvector.h"
  
  #include "libmesh/elem.h"
  
  using namespace libMesh;
  
  void assemble_stokes (EquationSystems& es,
                        const std::string& system_name);
  
  int main (int argc, char** argv)
  {
    LibMeshInit init (argc, argv);
  
    libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
  
    Mesh mesh(init.comm());
  
    MeshTools::Generation::build_square (mesh,
                                         15, 15,
                                         0., 1.,
                                         0., 1.,
                                         QUAD9);
  
    mesh.print_info();
  
    EquationSystems equation_systems (mesh);
  
    LinearImplicitSystem & system =
      equation_systems.add_system<LinearImplicitSystem> ("Stokes");
  
    system.add_variable ("u", SECOND);
    system.add_variable ("v", SECOND);
  
    system.add_variable ("p", FIRST);
  
    system.attach_assemble_function (assemble_stokes);
  
    equation_systems.init ();
  
    equation_systems.parameters.set<unsigned int>("linear solver maximum iterations") = 250;
    equation_systems.parameters.set<Real>        ("linear solver tolerance") = TOLERANCE;
  
    equation_systems.print_info();
  
    equation_systems.get_system("Stokes").solve();
  
  #ifdef LIBMESH_HAVE_EXODUS_API
    ExodusII_IO(mesh).write_equation_systems ("out.e",
                                        equation_systems);
  #endif // #ifdef LIBMESH_HAVE_EXODUS_API
  
    return 0;
  }
  
  void assemble_stokes (EquationSystems& es,
                        const std::string& system_name)
  {
    libmesh_assert_equal_to (system_name, "Stokes");
  
    const MeshBase& mesh = es.get_mesh();
  
    const unsigned int dim = mesh.mesh_dimension();
  
    LinearImplicitSystem & system =
      es.get_system<LinearImplicitSystem> ("Stokes");
  
    const unsigned int u_var = system.variable_number ("u");
    const unsigned int v_var = system.variable_number ("v");
    const unsigned int p_var = system.variable_number ("p");
  
    FEType fe_vel_type = system.variable_type(u_var);
  
    FEType fe_pres_type = system.variable_type(p_var);
  
    AutoPtr<FEBase> fe_vel  (FEBase::build(dim, fe_vel_type));
  
    AutoPtr<FEBase> fe_pres (FEBase::build(dim, fe_pres_type));
  
    QGauss qrule (dim, fe_vel_type.default_quadrature_order());
  
    fe_vel->attach_quadrature_rule (&qrule);
    fe_pres->attach_quadrature_rule (&qrule);
  
    const std::vector<Real>& JxW = fe_vel->get_JxW();
  
    const std::vector<std::vector<RealGradient> >& dphi = fe_vel->get_dphi();
  
    const std::vector<std::vector<Real> >& psi = fe_pres->get_phi();
  
    const DofMap & dof_map = system.get_dof_map();
  
    DenseMatrix<Number> Ke;
    DenseVector<Number> Fe;
  
    DenseSubMatrix<Number>
      Kuu(Ke), Kuv(Ke), Kup(Ke),
      Kvu(Ke), Kvv(Ke), Kvp(Ke),
      Kpu(Ke), Kpv(Ke), Kpp(Ke);
  
    DenseSubVector<Number>
      Fu(Fe),
      Fv(Fe),
      Fp(Fe);
  
    std::vector<dof_id_type> dof_indices;
    std::vector<dof_id_type> dof_indices_u;
    std::vector<dof_id_type> dof_indices_v;
    std::vector<dof_id_type> dof_indices_p;
  
  
    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);
        dof_map.dof_indices (elem, dof_indices_u, u_var);
        dof_map.dof_indices (elem, dof_indices_v, v_var);
        dof_map.dof_indices (elem, dof_indices_p, p_var);
  
        const unsigned int n_dofs   = dof_indices.size();
        const unsigned int n_u_dofs = dof_indices_u.size();
        const unsigned int n_v_dofs = dof_indices_v.size();
        const unsigned int n_p_dofs = dof_indices_p.size();
  
        fe_vel->reinit  (elem);
        fe_pres->reinit (elem);
  
        Ke.resize (n_dofs, n_dofs);
        Fe.resize (n_dofs);
  
        Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
        Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
        Kup.reposition (u_var*n_u_dofs, p_var*n_u_dofs, n_u_dofs, n_p_dofs);
  
        Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
        Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
        Kvp.reposition (v_var*n_v_dofs, p_var*n_v_dofs, n_v_dofs, n_p_dofs);
  
        Kpu.reposition (p_var*n_u_dofs, u_var*n_u_dofs, n_p_dofs, n_u_dofs);
        Kpv.reposition (p_var*n_u_dofs, v_var*n_u_dofs, n_p_dofs, n_v_dofs);
        Kpp.reposition (p_var*n_u_dofs, p_var*n_u_dofs, n_p_dofs, n_p_dofs);
  
        Fu.reposition (u_var*n_u_dofs, n_u_dofs);
        Fv.reposition (v_var*n_u_dofs, n_v_dofs);
        Fp.reposition (p_var*n_u_dofs, n_p_dofs);
  
        for (unsigned int qp=0; qp<qrule.n_points(); qp++)
          {
            for (unsigned int i=0; i<n_u_dofs; i++)
              for (unsigned int j=0; j<n_u_dofs; j++)
                Kuu(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
  
            for (unsigned int i=0; i<n_u_dofs; i++)
              for (unsigned int j=0; j<n_p_dofs; j++)
                Kup(i,j) += -JxW[qp]*psi[j][qp]*dphi[i][qp](0);
  
  
            for (unsigned int i=0; i<n_v_dofs; i++)
              for (unsigned int j=0; j<n_v_dofs; j++)
                Kvv(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
  
            for (unsigned int i=0; i<n_v_dofs; i++)
              for (unsigned int j=0; j<n_p_dofs; j++)
                Kvp(i,j) += -JxW[qp]*psi[j][qp]*dphi[i][qp](1);
  
  
            for (unsigned int i=0; i<n_p_dofs; i++)
              for (unsigned int j=0; j<n_u_dofs; j++)
                Kpu(i,j) += -JxW[qp]*psi[i][qp]*dphi[j][qp](0);
  
            for (unsigned int i=0; i<n_p_dofs; i++)
              for (unsigned int j=0; j<n_v_dofs; j++)
                Kpv(i,j) += -JxW[qp]*psi[i][qp]*dphi[j][qp](1);
  
          } // end of the quadrature point qp-loop
  
        {
          for (unsigned int s=0; s<elem->n_sides(); s++)
            if (elem->neighbor(s) == NULL)
              {
                AutoPtr<Elem> side (elem->build_side(s));
  
                for (unsigned int ns=0; ns<side->n_nodes(); ns++)
                  {
  
                    const Real yf = side->point(ns)(1);
  
                    const Real penalty = 1.e10;
  
  
                    const Real u_value = (yf > .99) ? 1. : 0.;
  
                    const Real v_value = 0.;
  
                    for (unsigned int n=0; n<elem->n_nodes(); n++)
                      if (elem->node(n) == side->node(ns))
                        {
                          Kuu(n,n) += penalty;
                          Kvv(n,n) += penalty;
  
                          Fu(n) += penalty*u_value;
                          Fv(n) += penalty*v_value;
                        }
                  } // end face node loop
              } // end if (elem->neighbor(side) == NULL)
        } // end boundary condition section
  
        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);
      } // end of element loop
  
    return;
  }



The console output of the program:

make[4]: Entering directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/systems_of_equations/systems_of_equations_ex1'
***************************************************************
* Running Example systems_of_equations_ex1:
*  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
***************************************************************
 
 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, "Stokes"
    Type "LinearImplicit"
    Variables={ "u" "v" } "p" 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" "LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" "CARTESIAN" 
    Approximation Orders="SECOND", "THIRD" "FIRST", "THIRD" 
    n_dofs()=2178
    n_local_dofs()=582
    n_constrained_dofs()=0
    n_local_constrained_dofs()=0
    n_vectors()=1
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 36.5546
      Average Off-Processor Bandwidth <= 2.70891
      Maximum  On-Processor Bandwidth <= 59
      Maximum Off-Processor Bandwidth <= 39
    DofMap Constraints
      Number of DoF Constraints = 0
      Number of Node Constraints = 0


 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:52:57 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.341842, Active time=0.239167                                                 |
 ----------------------------------------------------------------------------------------------------------------
| 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.0016      0.001614    0.0026      0.002571    0.67     1.07     |
|   build_sparsity()                 1         0.0014      0.001432    0.0035      0.003533    0.60     1.48     |
|   create_dof_constraints()         1         0.0002      0.000172    0.0002      0.000172    0.07     0.07     |
|   distribute_dofs()                1         0.0025      0.002501    0.0061      0.006056    1.05     2.53     |
|   dof_indices()                    479       0.0081      0.000017    0.0081      0.000017    3.40     3.40     |
|   prepare_send_list()              1         0.0000      0.000022    0.0000      0.000022    0.01     0.01     |
|   reinit()                         1         0.0032      0.003201    0.0032      0.003201    1.34     1.34     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          1         0.0006      0.000572    0.0028      0.002799    0.24     1.17     |
|                                                                                                                |
| ExodusII_IO                                                                                                    |
|   write_nodal_data()               1         0.1306      0.130605    0.1306      0.130605    54.61    54.61    |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        112       0.0005      0.000004    0.0005      0.000004    0.20     0.20     |
|   init_shape_functions()           2         0.0000      0.000023    0.0000      0.000023    0.02     0.02     |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             112       0.0006      0.000006    0.0006      0.000006    0.26     0.26     |
|   init_reference_to_physical_map() 2         0.0001      0.000047    0.0001      0.000047    0.04     0.04     |
|                                                                                                                |
| Mesh                                                                                                           |
|   find_neighbors()                 1         0.0007      0.000692    0.0008      0.000814    0.29     0.34     |
|   renumber_nodes_and_elem()        2         0.0001      0.000065    0.0001      0.000065    0.05     0.05     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   compute_hilbert_indices()        2         0.0020      0.000982    0.0020      0.000982    0.82     0.82     |
|   find_global_indices()            2         0.0003      0.000170    0.0030      0.001510    0.14     1.26     |
|   parallel_sort()                  2         0.0004      0.000176    0.0004      0.000217    0.15     0.18     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         1         0.0001      0.000072    0.1335      0.133549    0.03     55.84    |
|                                                                                                                |
| MeshTools::Generation                                                                                          |
|   build_cube()                     1         0.0005      0.000537    0.0005      0.000537    0.22     0.22     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      1         0.0026      0.002646    0.0041      0.004082    1.11     1.71     |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      9         0.0001      0.000011    0.0001      0.000014    0.04     0.05     |
|   max(bool)                        1         0.0000      0.000004    0.0000      0.000004    0.00     0.00     |
|   max(scalar)                      105       0.0005      0.000005    0.0005      0.000005    0.21     0.21     |
|   max(vector)                      24        0.0002      0.000007    0.0004      0.000019    0.07     0.19     |
|   min(bool)                        121       0.0005      0.000004    0.0005      0.000004    0.19     0.19     |
|   min(scalar)                      99        0.0007      0.000007    0.0007      0.000007    0.31     0.31     |
|   min(vector)                      24        0.0002      0.000010    0.0005      0.000023    0.10     0.23     |
|   probe()                          36        0.0002      0.000005    0.0002      0.000005    0.08     0.08     |
|   receive()                        36        0.0001      0.000004    0.0003      0.000009    0.05     0.13     |
|   send()                           36        0.0001      0.000003    0.0001      0.000003    0.04     0.04     |
|   send_receive()                   40        0.0002      0.000006    0.0007      0.000018    0.10     0.30     |
|   sum()                            20        0.0003      0.000014    0.0004      0.000018    0.12     0.15     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           36        0.0001      0.000002    0.0001      0.000002    0.02     0.02     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         1         0.0005      0.000467    0.0006      0.000608    0.20     0.25     |
|   set_parent_processor_ids()       1         0.0001      0.000070    0.0001      0.000070    0.03     0.03     |
|                                                                                                                |
| PetscLinearSolver                                                                                              |
|   solve()                          1         0.0728      0.072771    0.0728      0.072771    30.43    30.43    |
|                                                                                                                |
| System                                                                                                         |
|   assemble()                       1         0.0064      0.006443    0.0113      0.011332    2.69     4.74     |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            1318      0.2392                                          100.00            |
 ----------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example systems_of_equations_ex1:
*  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/systems_of_equations/systems_of_equations_ex1'

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

Hosted By:
SourceForge.net Logo