The source file systems_of_equations_ex6.C with comments:

Systems Example 6 - 3D Linear Elastic Cantilever

By David Knezevic

This is a 3D version of systems_of_equations_ex4. In this case we also compute and plot stresses.



C++ include files that we need
        #include <iostream>
        #include <algorithm>
        #include <math.h>
        
libMesh includes
        #include "libmesh/libmesh.h"
        #include "libmesh/mesh.h"
        #include "libmesh/mesh_generation.h"
        #include "libmesh/exodusII_io.h"
        #include "libmesh/gnuplot_io.h"
        #include "libmesh/linear_implicit_system.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_submatrix.h"
        #include "libmesh/dense_vector.h"
        #include "libmesh/dense_subvector.h"
        #include "libmesh/perf_log.h"
        #include "libmesh/elem.h"
        #include "libmesh/boundary_info.h"
        #include "libmesh/zero_function.h"
        #include "libmesh/dirichlet_boundaries.h"
        #include "libmesh/string_to_enum.h"
        #include "libmesh/getpot.h"
        
        #define x_scaling 1.3
        #define x_load 0.0
        #define y_load 0.0
        #define z_load -1.0
        
boundary IDs
        #define BOUNDARY_ID_MIN_Z 0
        #define BOUNDARY_ID_MIN_Y 1
        #define BOUNDARY_ID_MAX_X 2
        #define BOUNDARY_ID_MAX_Y 3
        #define BOUNDARY_ID_MIN_X 4
        #define BOUNDARY_ID_MAX_Z 5
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
Matrix and right-hand side assembly
        void assemble_elasticity(EquationSystems& es,
                                 const std::string& system_name);
        
Post-process the solution to compute stresses
        void compute_stresses(EquationSystems& es);
        
The Kronecker delta function, used in eval_elasticity_tensor
        Real kronecker_delta(unsigned int i,
                             unsigned int j);
        
Define the elasticity tensor, which is a fourth-order tensor i.e. it has four indices i,j,k,l
        Real eval_elasticity_tensor(unsigned int i,
                                    unsigned int j,
                                    unsigned int k,
                                    unsigned int l);
        
Begin the main program.
        int main (int argc, char** argv)
        {
Initialize libMesh and any dependent libaries
          LibMeshInit init (argc, argv);
        
Initialize the cantilever mesh
          const unsigned int dim = 3;
        
Make sure libMesh was compiled for 3D
          libmesh_example_assert(dim == LIBMESH_DIM, "3D support");
        
Create a 3D mesh distributed across the default MPI communicator.
          Mesh mesh(init.comm(), dim);
          MeshTools::Generation::build_cube (mesh,
                                             40,
                                             8,
                                             4,
                                             0., 1.*x_scaling,
                                             0., 0.3,
                                             0., 0.1,
                                             HEX8);
        
        
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 system named "Elasticity"
          LinearImplicitSystem& system =
            equation_systems.add_system<LinearImplicitSystem> ("Elasticity");
        
Add three displacement variables, u and v, to the system
          unsigned int u_var = system.add_variable("u", FIRST, LAGRANGE);
          unsigned int v_var = system.add_variable("v", FIRST, LAGRANGE);
          unsigned int w_var = system.add_variable("w", FIRST, LAGRANGE);
        
          system.attach_assemble_function (assemble_elasticity);
        
        
          std::set<boundary_id_type> boundary_ids;
          boundary_ids.insert(BOUNDARY_ID_MIN_X);
        
Create a vector storing the variable numbers which the BC applies to
          std::vector<unsigned int> variables;
          variables.push_back(u_var);
          variables.push_back(v_var);
          variables.push_back(w_var);
        
Create a ZeroFunction to initialize dirichlet_bc
          ZeroFunction<> zf;
        
          DirichletBoundary dirichlet_bc(boundary_ids,
                                         variables,
                                         &zf);
        
We must add the Dirichlet boundary condition _before_ we call equation_systems.init()
          system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
        
Also, initialize an ExplicitSystem to store stresses
          ExplicitSystem& stress_system =
            equation_systems.add_system<ExplicitSystem> ("StressSystem");
        
          stress_system.add_variable("sigma_00", CONSTANT, MONOMIAL);
          stress_system.add_variable("sigma_01", CONSTANT, MONOMIAL);
          stress_system.add_variable("sigma_02", CONSTANT, MONOMIAL);
          stress_system.add_variable("sigma_10", CONSTANT, MONOMIAL);
          stress_system.add_variable("sigma_11", CONSTANT, MONOMIAL);
          stress_system.add_variable("sigma_12", CONSTANT, MONOMIAL);
          stress_system.add_variable("sigma_20", CONSTANT, MONOMIAL);
          stress_system.add_variable("sigma_21", CONSTANT, MONOMIAL);
          stress_system.add_variable("sigma_22", CONSTANT, MONOMIAL);
          stress_system.add_variable("vonMises", CONSTANT, MONOMIAL);
        
Initialize the data structures for the equation system.
          equation_systems.init();
        
Print information about the system to the screen.
          equation_systems.print_info();
        
Solve the system
          system.solve();
        
Post-process the solution to compute the stresses
          compute_stresses(equation_systems);
        
Plot the solution
        #ifdef LIBMESH_HAVE_EXODUS_API
          ExodusII_IO (mesh).write_equation_systems("displacement.e",equation_systems);
        #endif // #ifdef LIBMESH_HAVE_EXODUS_API
        
All done.
          return 0;
        }
        
        
        void assemble_elasticity(EquationSystems& es,
                                 const std::string& system_name)
        {
          libmesh_assert_equal_to (system_name, "Elasticity");
        
          const MeshBase& mesh = es.get_mesh();
        
          const unsigned int dim = mesh.mesh_dimension();
        
          LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Elasticity");
        
          const unsigned int n_components = 3;
          const unsigned int u_var = system.variable_number ("u");
          const unsigned int v_var = system.variable_number ("v");
          const unsigned int w_var = system.variable_number ("w");
        
          const DofMap& dof_map = system.get_dof_map();
          FEType fe_type = dof_map.variable_type(u_var);
          AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
          QGauss qrule (dim, fe_type.default_quadrature_order());
          fe->attach_quadrature_rule (&qrule);
        
          AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
          QGauss qface(dim-1, fe_type.default_quadrature_order());
          fe_face->attach_quadrature_rule (&qface);
        
          const std::vector<Real>& JxW = fe->get_JxW();
          const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
        
          DenseMatrix<Number> Ke;
          DenseVector<Number> Fe;
        
          DenseSubMatrix<Number>
            Kuu(Ke), Kuv(Ke), Kuw(Ke),
            Kvu(Ke), Kvv(Ke), Kvw(Ke),
            Kwu(Ke), Kwv(Ke), Kww(Ke);
        
          DenseSubVector<Number>
            Fu(Fe),
            Fv(Fe),
            Fw(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_w;
        
          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_w, w_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_w_dofs = dof_indices_w.size();
        
              fe->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);
              Kuw.reposition (u_var*n_u_dofs, w_var*n_u_dofs, n_u_dofs, n_w_dofs);
        
              Kvu.reposition (v_var*n_u_dofs, u_var*n_u_dofs, n_v_dofs, n_u_dofs);
              Kvv.reposition (v_var*n_u_dofs, v_var*n_u_dofs, n_v_dofs, n_v_dofs);
              Kvw.reposition (v_var*n_u_dofs, w_var*n_u_dofs, n_v_dofs, n_w_dofs);
        
              Kwu.reposition (w_var*n_u_dofs, u_var*n_u_dofs, n_w_dofs, n_u_dofs);
              Kwv.reposition (w_var*n_u_dofs, v_var*n_u_dofs, n_w_dofs, n_v_dofs);
              Kww.reposition (w_var*n_u_dofs, w_var*n_u_dofs, n_w_dofs, n_w_dofs);
        
              Fu.reposition (u_var*n_u_dofs, n_u_dofs);
              Fv.reposition (v_var*n_u_dofs, n_v_dofs);
              Fw.reposition (w_var*n_u_dofs, n_w_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++)
                    {
Tensor indices
                      unsigned int C_i=0, C_k=0;
        
                      for(unsigned int C_j=0; C_j<n_components; C_j++)
                        for(unsigned int C_l=0; C_l<n_components; C_l++)
                        {
                          Kuu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                        }
                    }
        
                  for (unsigned int i=0; i<n_u_dofs; i++)
                    for (unsigned int j=0; j<n_v_dofs; j++)
                    {
Tensor indices
                      unsigned int C_i=0, C_k=1;
        
                      for(unsigned int C_j=0; C_j<n_components; C_j++)
                        for(unsigned int C_l=0; C_l<n_components; C_l++)
                        {
                          Kuv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                        }
                    }
        
                  for (unsigned int i=0; i<n_u_dofs; i++)
                    for (unsigned int j=0; j<n_w_dofs; j++)
                    {
Tensor indices
                      unsigned int C_i=0, C_k=2;
        
                      for(unsigned int C_j=0; C_j<n_components; C_j++)
                        for(unsigned int C_l=0; C_l<n_components; C_l++)
                        {
                          Kuw(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                        }
                    }
        
                  for (unsigned int i=0; i<n_v_dofs; i++)
                    for (unsigned int j=0; j<n_u_dofs; j++)
                    {
Tensor indices
                      unsigned int C_i = 1, C_k = 0;
        
                      for(unsigned int C_j=0; C_j<n_components; C_j++)
                        for(unsigned int C_l=0; C_l<n_components; C_l++)
                        {
                          Kvu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                        }
                    }
        
                  for (unsigned int i=0; i<n_v_dofs; i++)
                    for (unsigned int j=0; j<n_v_dofs; j++)
                    {
Tensor indices
                      unsigned int C_i = 1, C_k = 1;
        
                      for(unsigned int C_j=0; C_j<n_components; C_j++)
                        for(unsigned int C_l=0; C_l<n_components; C_l++)
                        {
                          Kvv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                        }
                    }
        
                  for (unsigned int i=0; i<n_v_dofs; i++)
                    for (unsigned int j=0; j<n_w_dofs; j++)
                    {
Tensor indices
                      unsigned int C_i = 1, C_k = 2;
        
                      for(unsigned int C_j=0; C_j<n_components; C_j++)
                        for(unsigned int C_l=0; C_l<n_components; C_l++)
                        {
                          Kvw(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                        }
                    }
        
                  for (unsigned int i=0; i<n_w_dofs; i++)
                    for (unsigned int j=0; j<n_u_dofs; j++)
                    {
Tensor indices
                      unsigned int C_i = 2, C_k = 0;
        
                      for(unsigned int C_j=0; C_j<n_components; C_j++)
                        for(unsigned int C_l=0; C_l<n_components; C_l++)
                        {
                          Kwu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                        }
                    }
        
                  for (unsigned int i=0; i<n_w_dofs; i++)
                    for (unsigned int j=0; j<n_v_dofs; j++)
                    {
Tensor indices
                      unsigned int C_i = 2, C_k = 1;
        
                      for(unsigned int C_j=0; C_j<n_components; C_j++)
                        for(unsigned int C_l=0; C_l<n_components; C_l++)
                        {
                          Kwv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                        }
                    }
        
                  for (unsigned int i=0; i<n_w_dofs; i++)
                    for (unsigned int j=0; j<n_w_dofs; j++)
                    {
Tensor indices
                      unsigned int C_i = 2, C_k = 2;
        
                      for(unsigned int C_j=0; C_j<n_components; C_j++)
                        for(unsigned int C_l=0; C_l<n_components; C_l++)
                        {
                          Kww(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                        }
                    }
              }
        
              {
                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++)
                      {
Apply a traction
                        if( mesh.boundary_info->has_boundary_id(elem, side, BOUNDARY_ID_MAX_X) )
                        {
                          for (unsigned int i=0; i<n_v_dofs; i++)
                          {
                            Fu(i) += JxW_face[qp] * x_load * phi_face[i][qp];
                          }
                          for (unsigned int i=0; i<n_v_dofs; i++)
                          {
                            Fv(i) += JxW_face[qp] * y_load * phi_face[i][qp];
                          }
                          for (unsigned int i=0; i<n_v_dofs; i++)
                          {
                            Fw(i) += JxW_face[qp] * z_load * 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);
            }
        }
        
        void compute_stresses(EquationSystems& es)
        {
          const MeshBase& mesh = es.get_mesh();
        
          const unsigned int dim = mesh.mesh_dimension();
        
          LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Elasticity");
        
          unsigned int displacement_vars[3];
          displacement_vars[0] = system.variable_number ("u");
          displacement_vars[1] = system.variable_number ("v");
          displacement_vars[2] = system.variable_number ("w");
          const unsigned int u_var = system.variable_number ("u");
        
          const DofMap& dof_map = system.get_dof_map();
          FEType fe_type = dof_map.variable_type(u_var);
          AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
          QGauss qrule (dim, fe_type.default_quadrature_order());
          fe->attach_quadrature_rule (&qrule);
        
          const std::vector<Real>& JxW = fe->get_JxW();
          const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
        
Also, get a reference to the ExplicitSystem
          ExplicitSystem& stress_system = es.get_system<ExplicitSystem>("StressSystem");
          const DofMap& stress_dof_map = stress_system.get_dof_map();
          unsigned int sigma_vars[3][3];
          sigma_vars[0][0] = stress_system.variable_number ("sigma_00");
          sigma_vars[0][1] = stress_system.variable_number ("sigma_01");
          sigma_vars[0][2] = stress_system.variable_number ("sigma_02");
          sigma_vars[1][0] = stress_system.variable_number ("sigma_10");
          sigma_vars[1][1] = stress_system.variable_number ("sigma_11");
          sigma_vars[1][2] = stress_system.variable_number ("sigma_12");
          sigma_vars[2][0] = stress_system.variable_number ("sigma_20");
          sigma_vars[2][1] = stress_system.variable_number ("sigma_21");
          sigma_vars[2][2] = stress_system.variable_number ("sigma_22");
          unsigned int vonMises_var = stress_system.variable_number ("vonMises");
        
Storage for the stress dof indices on each element
          std::vector< std::vector<dof_id_type> > dof_indices_var(system.n_vars());
          std::vector<dof_id_type> stress_dof_indices_var;
        
To store the stress tensor on each element
          DenseMatrix<Number> elem_sigma;
        
          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;
        
            for(unsigned int var=0; var<3; var++)
            {
              dof_map.dof_indices (elem, dof_indices_var[var], displacement_vars[var]);
            }
        
            fe->reinit (elem);
        
            elem_sigma.resize(3,3);
        
            for (unsigned int qp=0; qp<qrule.n_points(); qp++)
            {
              for(unsigned int C_i=0; C_i<3; C_i++)
                for(unsigned int C_j=0; C_j<3; C_j++)
                  for(unsigned int C_k=0; C_k<3; C_k++)
                  {
                    const unsigned int n_var_dofs = dof_indices_var[C_k].size();
        
Get the gradient at this quadrature point
                    Gradient displacement_gradient;
                    for(unsigned int l=0; l<n_var_dofs; l++)
                    {
                      displacement_gradient.add_scaled(dphi[l][qp], system.current_solution(dof_indices_var[C_k][l]));
                    }
        
                    for(unsigned int C_l=0; C_l<3; C_l++)
                    {
                      elem_sigma(C_i,C_j) += JxW[qp]*( eval_elasticity_tensor(C_i,C_j,C_k,C_l) * displacement_gradient(C_l) );
                    }
        
                  }
            }
        
Get the average stresses by dividing by the element volume
            elem_sigma.scale(1./elem->volume());
        
load elem_sigma data into stress_system
            for(unsigned int i=0; i<3; i++)
              for(unsigned int j=0; j<3; j++)
              {
                stress_dof_map.dof_indices (elem, stress_dof_indices_var, sigma_vars[i][j]);
        
We are using CONSTANT MONOMIAL basis functions, hence we only need to get one dof index per variable
                dof_id_type dof_index = stress_dof_indices_var[0];
        
                if( (stress_system.solution->first_local_index() <= dof_index) &&
                    (dof_index < stress_system.solution->last_local_index()) )
                {
                  stress_system.solution->set(dof_index, elem_sigma(i,j));
                }
        
              }
        
Also, the von Mises stress
            Number vonMises_value = std::sqrt( 0.5*( pow(elem_sigma(0,0) - elem_sigma(1,1),2.) +
                                                     pow(elem_sigma(1,1) - elem_sigma(2,2),2.) +
                                                     pow(elem_sigma(2,2) - elem_sigma(0,0),2.) +
                                                     6.*(pow(elem_sigma(0,1),2.) + pow(elem_sigma(1,2),2.) + pow(elem_sigma(2,0),2.))
                                                   ) );
            stress_dof_map.dof_indices (elem, stress_dof_indices_var, vonMises_var);
            dof_id_type dof_index = stress_dof_indices_var[0];
            if( (stress_system.solution->first_local_index() <= dof_index) &&
                (dof_index < stress_system.solution->last_local_index()) )
            {
              stress_system.solution->set(dof_index, vonMises_value);
            }
        
          }
        
Should call close and update when we set vector entries directly
          stress_system.solution->close();
          stress_system.update();
        }
        
        Real kronecker_delta(unsigned int i,
                             unsigned int j)
        {
          return i == j ? 1. : 0.;
        }
        
        Real eval_elasticity_tensor(unsigned int i,
                                    unsigned int j,
                                    unsigned int k,
                                    unsigned int l)
        {
Define the Poisson ratio and Young's modulus
          const Real nu = 0.3;
          const Real E  = 1.;
        
Define the Lame constants (lambda_1 and lambda_2) based on nu and E
          const Real lambda_1 = E * nu / ( (1. + nu) * (1. - 2.*nu) );
          const Real lambda_2 = 0.5 * E / (1. + nu);
        
          return lambda_1 * kronecker_delta(i,j) * kronecker_delta(k,l)
               + lambda_2 * (kronecker_delta(i,k) * kronecker_delta(j,l) + kronecker_delta(i,l) * kronecker_delta(j,k));
        }



The source file systems_of_equations_ex6.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/gnuplot_io.h"
  #include "libmesh/linear_implicit_system.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_submatrix.h"
  #include "libmesh/dense_vector.h"
  #include "libmesh/dense_subvector.h"
  #include "libmesh/perf_log.h"
  #include "libmesh/elem.h"
  #include "libmesh/boundary_info.h"
  #include "libmesh/zero_function.h"
  #include "libmesh/dirichlet_boundaries.h"
  #include "libmesh/string_to_enum.h"
  #include "libmesh/getpot.h"
  
  #define x_scaling 1.3
  #define x_load 0.0
  #define y_load 0.0
  #define z_load -1.0
  
  #define BOUNDARY_ID_MIN_Z 0
  #define BOUNDARY_ID_MIN_Y 1
  #define BOUNDARY_ID_MAX_X 2
  #define BOUNDARY_ID_MAX_Y 3
  #define BOUNDARY_ID_MIN_X 4
  #define BOUNDARY_ID_MAX_Z 5
  
  using namespace libMesh;
  
  void assemble_elasticity(EquationSystems& es,
                           const std::string& system_name);
  
  void compute_stresses(EquationSystems& es);
  
  Real kronecker_delta(unsigned int i,
                       unsigned int j);
  
  Real eval_elasticity_tensor(unsigned int i,
                              unsigned int j,
                              unsigned int k,
                              unsigned int l);
  
  int main (int argc, char** argv)
  {
    LibMeshInit init (argc, argv);
  
    const unsigned int dim = 3;
  
    libmesh_example_assert(dim == LIBMESH_DIM, "3D support");
  
    Mesh mesh(init.comm(), dim);
    MeshTools::Generation::build_cube (mesh,
                                       40,
                                       8,
                                       4,
                                       0., 1.*x_scaling,
                                       0., 0.3,
                                       0., 0.1,
                                       HEX8);
  
  
    mesh.print_info();
  
  
    EquationSystems equation_systems (mesh);
  
    LinearImplicitSystem& system =
      equation_systems.add_system<LinearImplicitSystem> ("Elasticity");
  
    unsigned int u_var = system.add_variable("u", FIRST, LAGRANGE);
    unsigned int v_var = system.add_variable("v", FIRST, LAGRANGE);
    unsigned int w_var = system.add_variable("w", FIRST, LAGRANGE);
  
    system.attach_assemble_function (assemble_elasticity);
  
  
    std::set<boundary_id_type> boundary_ids;
    boundary_ids.insert(BOUNDARY_ID_MIN_X);
  
    std::vector<unsigned int> variables;
    variables.push_back(u_var);
    variables.push_back(v_var);
    variables.push_back(w_var);
  
    ZeroFunction<> zf;
  
    DirichletBoundary dirichlet_bc(boundary_ids,
                                   variables,
                                   &zf);
  
    system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
  
    ExplicitSystem& stress_system =
      equation_systems.add_system<ExplicitSystem> ("StressSystem");
  
    stress_system.add_variable("sigma_00", CONSTANT, MONOMIAL);
    stress_system.add_variable("sigma_01", CONSTANT, MONOMIAL);
    stress_system.add_variable("sigma_02", CONSTANT, MONOMIAL);
    stress_system.add_variable("sigma_10", CONSTANT, MONOMIAL);
    stress_system.add_variable("sigma_11", CONSTANT, MONOMIAL);
    stress_system.add_variable("sigma_12", CONSTANT, MONOMIAL);
    stress_system.add_variable("sigma_20", CONSTANT, MONOMIAL);
    stress_system.add_variable("sigma_21", CONSTANT, MONOMIAL);
    stress_system.add_variable("sigma_22", CONSTANT, MONOMIAL);
    stress_system.add_variable("vonMises", CONSTANT, MONOMIAL);
  
    equation_systems.init();
  
    equation_systems.print_info();
  
    system.solve();
  
    compute_stresses(equation_systems);
  
  #ifdef LIBMESH_HAVE_EXODUS_API
    ExodusII_IO (mesh).write_equation_systems("displacement.e",equation_systems);
  #endif // #ifdef LIBMESH_HAVE_EXODUS_API
  
    return 0;
  }
  
  
  void assemble_elasticity(EquationSystems& es,
                           const std::string& system_name)
  {
    libmesh_assert_equal_to (system_name, "Elasticity");
  
    const MeshBase& mesh = es.get_mesh();
  
    const unsigned int dim = mesh.mesh_dimension();
  
    LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Elasticity");
  
    const unsigned int n_components = 3;
    const unsigned int u_var = system.variable_number ("u");
    const unsigned int v_var = system.variable_number ("v");
    const unsigned int w_var = system.variable_number ("w");
  
    const DofMap& dof_map = system.get_dof_map();
    FEType fe_type = dof_map.variable_type(u_var);
    AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
    QGauss qrule (dim, fe_type.default_quadrature_order());
    fe->attach_quadrature_rule (&qrule);
  
    AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
    QGauss qface(dim-1, fe_type.default_quadrature_order());
    fe_face->attach_quadrature_rule (&qface);
  
    const std::vector<Real>& JxW = fe->get_JxW();
    const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
  
    DenseMatrix<Number> Ke;
    DenseVector<Number> Fe;
  
    DenseSubMatrix<Number>
      Kuu(Ke), Kuv(Ke), Kuw(Ke),
      Kvu(Ke), Kvv(Ke), Kvw(Ke),
      Kwu(Ke), Kwv(Ke), Kww(Ke);
  
    DenseSubVector<Number>
      Fu(Fe),
      Fv(Fe),
      Fw(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_w;
  
    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_w, w_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_w_dofs = dof_indices_w.size();
  
        fe->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);
        Kuw.reposition (u_var*n_u_dofs, w_var*n_u_dofs, n_u_dofs, n_w_dofs);
  
        Kvu.reposition (v_var*n_u_dofs, u_var*n_u_dofs, n_v_dofs, n_u_dofs);
        Kvv.reposition (v_var*n_u_dofs, v_var*n_u_dofs, n_v_dofs, n_v_dofs);
        Kvw.reposition (v_var*n_u_dofs, w_var*n_u_dofs, n_v_dofs, n_w_dofs);
  
        Kwu.reposition (w_var*n_u_dofs, u_var*n_u_dofs, n_w_dofs, n_u_dofs);
        Kwv.reposition (w_var*n_u_dofs, v_var*n_u_dofs, n_w_dofs, n_v_dofs);
        Kww.reposition (w_var*n_u_dofs, w_var*n_u_dofs, n_w_dofs, n_w_dofs);
  
        Fu.reposition (u_var*n_u_dofs, n_u_dofs);
        Fv.reposition (v_var*n_u_dofs, n_v_dofs);
        Fw.reposition (w_var*n_u_dofs, n_w_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++)
              {
                unsigned int C_i=0, C_k=0;
  
                for(unsigned int C_j=0; C_j<n_components; C_j++)
                  for(unsigned int C_l=0; C_l<n_components; C_l++)
                  {
                    Kuu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                  }
              }
  
            for (unsigned int i=0; i<n_u_dofs; i++)
              for (unsigned int j=0; j<n_v_dofs; j++)
              {
                unsigned int C_i=0, C_k=1;
  
                for(unsigned int C_j=0; C_j<n_components; C_j++)
                  for(unsigned int C_l=0; C_l<n_components; C_l++)
                  {
                    Kuv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                  }
              }
  
            for (unsigned int i=0; i<n_u_dofs; i++)
              for (unsigned int j=0; j<n_w_dofs; j++)
              {
                unsigned int C_i=0, C_k=2;
  
                for(unsigned int C_j=0; C_j<n_components; C_j++)
                  for(unsigned int C_l=0; C_l<n_components; C_l++)
                  {
                    Kuw(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                  }
              }
  
            for (unsigned int i=0; i<n_v_dofs; i++)
              for (unsigned int j=0; j<n_u_dofs; j++)
              {
                unsigned int C_i = 1, C_k = 0;
  
                for(unsigned int C_j=0; C_j<n_components; C_j++)
                  for(unsigned int C_l=0; C_l<n_components; C_l++)
                  {
                    Kvu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                  }
              }
  
            for (unsigned int i=0; i<n_v_dofs; i++)
              for (unsigned int j=0; j<n_v_dofs; j++)
              {
                unsigned int C_i = 1, C_k = 1;
  
                for(unsigned int C_j=0; C_j<n_components; C_j++)
                  for(unsigned int C_l=0; C_l<n_components; C_l++)
                  {
                    Kvv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                  }
              }
  
            for (unsigned int i=0; i<n_v_dofs; i++)
              for (unsigned int j=0; j<n_w_dofs; j++)
              {
                unsigned int C_i = 1, C_k = 2;
  
                for(unsigned int C_j=0; C_j<n_components; C_j++)
                  for(unsigned int C_l=0; C_l<n_components; C_l++)
                  {
                    Kvw(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                  }
              }
  
            for (unsigned int i=0; i<n_w_dofs; i++)
              for (unsigned int j=0; j<n_u_dofs; j++)
              {
                unsigned int C_i = 2, C_k = 0;
  
                for(unsigned int C_j=0; C_j<n_components; C_j++)
                  for(unsigned int C_l=0; C_l<n_components; C_l++)
                  {
                    Kwu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                  }
              }
  
            for (unsigned int i=0; i<n_w_dofs; i++)
              for (unsigned int j=0; j<n_v_dofs; j++)
              {
                unsigned int C_i = 2, C_k = 1;
  
                for(unsigned int C_j=0; C_j<n_components; C_j++)
                  for(unsigned int C_l=0; C_l<n_components; C_l++)
                  {
                    Kwv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                  }
              }
  
            for (unsigned int i=0; i<n_w_dofs; i++)
              for (unsigned int j=0; j<n_w_dofs; j++)
              {
                unsigned int C_i = 2, C_k = 2;
  
                for(unsigned int C_j=0; C_j<n_components; C_j++)
                  for(unsigned int C_l=0; C_l<n_components; C_l++)
                  {
                    Kww(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                  }
              }
        }
  
        {
          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++)
                {
                  if( mesh.boundary_info->has_boundary_id(elem, side, BOUNDARY_ID_MAX_X) )
                  {
                    for (unsigned int i=0; i<n_v_dofs; i++)
                    {
                      Fu(i) += JxW_face[qp] * x_load * phi_face[i][qp];
                    }
                    for (unsigned int i=0; i<n_v_dofs; i++)
                    {
                      Fv(i) += JxW_face[qp] * y_load * phi_face[i][qp];
                    }
                    for (unsigned int i=0; i<n_v_dofs; i++)
                    {
                      Fw(i) += JxW_face[qp] * z_load * 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);
      }
  }
  
  void compute_stresses(EquationSystems& es)
  {
    const MeshBase& mesh = es.get_mesh();
  
    const unsigned int dim = mesh.mesh_dimension();
  
    LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Elasticity");
  
    unsigned int displacement_vars[3];
    displacement_vars[0] = system.variable_number ("u");
    displacement_vars[1] = system.variable_number ("v");
    displacement_vars[2] = system.variable_number ("w");
    const unsigned int u_var = system.variable_number ("u");
  
    const DofMap& dof_map = system.get_dof_map();
    FEType fe_type = dof_map.variable_type(u_var);
    AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
    QGauss qrule (dim, fe_type.default_quadrature_order());
    fe->attach_quadrature_rule (&qrule);
  
    const std::vector<Real>& JxW = fe->get_JxW();
    const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
  
    ExplicitSystem& stress_system = es.get_system<ExplicitSystem>("StressSystem");
    const DofMap& stress_dof_map = stress_system.get_dof_map();
    unsigned int sigma_vars[3][3];
    sigma_vars[0][0] = stress_system.variable_number ("sigma_00");
    sigma_vars[0][1] = stress_system.variable_number ("sigma_01");
    sigma_vars[0][2] = stress_system.variable_number ("sigma_02");
    sigma_vars[1][0] = stress_system.variable_number ("sigma_10");
    sigma_vars[1][1] = stress_system.variable_number ("sigma_11");
    sigma_vars[1][2] = stress_system.variable_number ("sigma_12");
    sigma_vars[2][0] = stress_system.variable_number ("sigma_20");
    sigma_vars[2][1] = stress_system.variable_number ("sigma_21");
    sigma_vars[2][2] = stress_system.variable_number ("sigma_22");
    unsigned int vonMises_var = stress_system.variable_number ("vonMises");
  
    std::vector< std::vector<dof_id_type> > dof_indices_var(system.n_vars());
    std::vector<dof_id_type> stress_dof_indices_var;
  
    DenseMatrix<Number> elem_sigma;
  
    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;
  
      for(unsigned int var=0; var<3; var++)
      {
        dof_map.dof_indices (elem, dof_indices_var[var], displacement_vars[var]);
      }
  
      fe->reinit (elem);
  
      elem_sigma.resize(3,3);
  
      for (unsigned int qp=0; qp<qrule.n_points(); qp++)
      {
        for(unsigned int C_i=0; C_i<3; C_i++)
          for(unsigned int C_j=0; C_j<3; C_j++)
            for(unsigned int C_k=0; C_k<3; C_k++)
            {
              const unsigned int n_var_dofs = dof_indices_var[C_k].size();
  
              Gradient displacement_gradient;
              for(unsigned int l=0; l<n_var_dofs; l++)
              {
                displacement_gradient.add_scaled(dphi[l][qp], system.current_solution(dof_indices_var[C_k][l]));
              }
  
              for(unsigned int C_l=0; C_l<3; C_l++)
              {
                elem_sigma(C_i,C_j) += JxW[qp]*( eval_elasticity_tensor(C_i,C_j,C_k,C_l) * displacement_gradient(C_l) );
              }
  
            }
      }
  
      elem_sigma.scale(1./elem->volume());
  
      for(unsigned int i=0; i<3; i++)
        for(unsigned int j=0; j<3; j++)
        {
          stress_dof_map.dof_indices (elem, stress_dof_indices_var, sigma_vars[i][j]);
  
          dof_id_type dof_index = stress_dof_indices_var[0];
  
          if( (stress_system.solution->first_local_index() <= dof_index) &&
              (dof_index < stress_system.solution->last_local_index()) )
          {
            stress_system.solution->set(dof_index, elem_sigma(i,j));
          }
  
        }
  
      Number vonMises_value = std::sqrt( 0.5*( pow(elem_sigma(0,0) - elem_sigma(1,1),2.) +
                                               pow(elem_sigma(1,1) - elem_sigma(2,2),2.) +
                                               pow(elem_sigma(2,2) - elem_sigma(0,0),2.) +
                                               6.*(pow(elem_sigma(0,1),2.) + pow(elem_sigma(1,2),2.) + pow(elem_sigma(2,0),2.))
                                             ) );
      stress_dof_map.dof_indices (elem, stress_dof_indices_var, vonMises_var);
      dof_id_type dof_index = stress_dof_indices_var[0];
      if( (stress_system.solution->first_local_index() <= dof_index) &&
          (dof_index < stress_system.solution->last_local_index()) )
      {
        stress_system.solution->set(dof_index, vonMises_value);
      }
  
    }
  
    stress_system.solution->close();
    stress_system.update();
  }
  
  Real kronecker_delta(unsigned int i,
                       unsigned int j)
  {
    return i == j ? 1. : 0.;
  }
  
  Real eval_elasticity_tensor(unsigned int i,
                              unsigned int j,
                              unsigned int k,
                              unsigned int l)
  {
    const Real nu = 0.3;
    const Real E  = 1.;
  
    const Real lambda_1 = E * nu / ( (1. + nu) * (1. - 2.*nu) );
    const Real lambda_2 = 0.5 * E / (1. + nu);
  
    return lambda_1 * kronecker_delta(i,j) * kronecker_delta(k,l)
         + lambda_2 * (kronecker_delta(i,k) * kronecker_delta(j,l) + kronecker_delta(i,l) * kronecker_delta(j,k));
  }



The console output of the program:

make[4]: Entering directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/systems_of_equations/systems_of_equations_ex6'
***************************************************************
* Running Example systems_of_equations_ex6:
*  mpirun -np 4 example-devel -ksp_type cg -pc_type bjacobi -sub_pc_type icc -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
 
 Mesh Information:
  mesh_dimension()=3
  spatial_dimension()=3
  n_nodes()=1845
    n_local_nodes()=495
  n_elem()=1280
    n_local_elem()=320
    n_active_elem()=1280
  n_subdomains()=1
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0

 EquationSystems
  n_systems()=2
   System #0, "Elasticity"
    Type "LinearImplicit"
    Variables={ "u" "v" "w" } 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="FIRST", "THIRD" 
    n_dofs()=5535
    n_local_dofs()=1485
    n_constrained_dofs()=135
    n_local_constrained_dofs()=135
    n_vectors()=1
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 62.3577
      Average Off-Processor Bandwidth <= 3.17073
      Maximum  On-Processor Bandwidth <= 81
      Maximum Off-Processor Bandwidth <= 27
    DofMap Constraints
      Number of DoF Constraints = 135
      Average DoF Constraint Length= 0
      Number of Node Constraints = 0
   System #1, "StressSystem"
    Type "Explicit"
    Variables={ "sigma_00" "sigma_01" "sigma_02" "sigma_10" "sigma_11" "sigma_12" "sigma_20" "sigma_21" "sigma_22" "vonMises" } 
    Finite Element Types="MONOMIAL", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="CONSTANT", "THIRD" 
    n_dofs()=12800
    n_local_dofs()=3200
    n_constrained_dofs()=0
    n_local_constrained_dofs()=0
    n_vectors()=1
    n_matrices()=0
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 0
      Average Off-Processor Bandwidth <= 0
      Maximum  On-Processor Bandwidth <= 0
      Maximum Off-Processor Bandwidth <= 0
    DofMap Constraints
      Number of DoF Constraints = 0
      Number of Node Constraints = 0


 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:54:27 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=1.45675, Active time=1.37421                                                   |
 ----------------------------------------------------------------------------------------------------------------
| 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()     2         0.0073      0.003625    0.0087      0.004344    0.53     0.63     |
|   build_constraint_matrix()        320       0.0007      0.000002    0.0007      0.000002    0.05     0.05     |
|   build_sparsity()                 1         0.0059      0.005939    0.0159      0.015914    0.43     1.16     |
|   cnstrn_elem_mat_vec()            320       0.0047      0.000015    0.0047      0.000015    0.34     0.34     |
|   create_dof_constraints()         2         0.0080      0.003997    0.0258      0.012877    0.58     1.87     |
|   distribute_dofs()                2         0.0043      0.002128    0.0209      0.010448    0.31     1.52     |
|   dof_indices()                    13888     0.0495      0.000004    0.0495      0.000004    3.60     3.60     |
|   prepare_send_list()              2         0.0001      0.000033    0.0001      0.000033    0.00     0.00     |
|   reinit()                         2         0.0080      0.004017    0.0080      0.004017    0.58     0.58     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          1         0.0053      0.005337    0.0279      0.027888    0.39     2.03     |
|                                                                                                                |
| ExodusII_IO                                                                                                    |
|   write_nodal_data()               1         0.3879      0.387852    0.3879      0.387852    28.22    28.22    |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        912       0.0027      0.000003    0.0027      0.000003    0.20     0.20     |
|   init_shape_functions()           274       0.0004      0.000001    0.0004      0.000001    0.03     0.03     |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             912       0.0032      0.000004    0.0032      0.000004    0.23     0.23     |
|   compute_face_map()               272       0.0010      0.000004    0.0010      0.000004    0.07     0.07     |
|   init_face_shape_functions()      1         0.0000      0.000014    0.0000      0.000014    0.00     0.00     |
|   init_reference_to_physical_map() 274       0.0027      0.000010    0.0027      0.000010    0.19     0.19     |
|                                                                                                                |
| Mesh                                                                                                           |
|   find_neighbors()                 1         0.0034      0.003359    0.0058      0.005814    0.24     0.42     |
|   renumber_nodes_and_elem()        2         0.0003      0.000138    0.0003      0.000138    0.02     0.02     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   compute_hilbert_indices()        2         0.0062      0.003075    0.0062      0.003075    0.45     0.45     |
|   find_global_indices()            2         0.0010      0.000503    0.0130      0.006517    0.07     0.95     |
|   parallel_sort()                  2         0.0005      0.000230    0.0054      0.002700    0.03     0.39     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         1         0.0001      0.000074    0.4159      0.415894    0.01     30.26    |
|                                                                                                                |
| MeshTools::Generation                                                                                          |
|   build_cube()                     1         0.0011      0.001078    0.0011      0.001078    0.08     0.08     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      1         0.0102      0.010169    0.0170      0.016996    0.74     1.24     |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      12        0.0060      0.000498    0.0061      0.000509    0.44     0.44     |
|   max(bool)                        1         0.0000      0.000003    0.0000      0.000003    0.00     0.00     |
|   max(scalar)                      137       0.0007      0.000005    0.0007      0.000005    0.05     0.05     |
|   max(vector)                      31        0.0002      0.000005    0.0005      0.000018    0.01     0.04     |
|   min(bool)                        157       0.0006      0.000004    0.0006      0.000004    0.04     0.04     |
|   min(scalar)                      128       0.0557      0.000435    0.0557      0.000435    4.05     4.05     |
|   min(vector)                      31        0.0002      0.000008    0.0010      0.000033    0.02     0.07     |
|   probe()                          48        0.0019      0.000040    0.0019      0.000040    0.14     0.14     |
|   receive()                        48        0.0002      0.000003    0.0021      0.000044    0.01     0.15     |
|   send()                           48        0.0001      0.000002    0.0001      0.000002    0.01     0.01     |
|   send_receive()                   52        0.0003      0.000005    0.0025      0.000048    0.02     0.18     |
|   sum()                            28        0.0063      0.000225    0.0135      0.000481    0.46     0.98     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           48        0.0001      0.000001    0.0001      0.000001    0.00     0.00     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         1         0.0008      0.000815    0.0065      0.006537    0.06     0.48     |
|   set_parent_processor_ids()       1         0.0002      0.000205    0.0002      0.000205    0.01     0.01     |
|                                                                                                                |
| PetscLinearSolver                                                                                              |
|   solve()                          1         0.5068      0.506793    0.5068      0.506793    36.88    36.88    |
|                                                                                                                |
| System                                                                                                         |
|   assemble()                       1         0.2800      0.279984    0.3021      0.302119    20.37    21.98    |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            17971     1.3742                                          100.00            |
 ----------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example systems_of_equations_ex6:
*  mpirun -np 4 example-devel -ksp_type cg -pc_type bjacobi -sub_pc_type icc -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_ex6'

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

Hosted By:
SourceForge.net Logo