The source file miscellaneous_ex4.C with comments:

Miscellaneous Example 4 - Using a shell matrix



This example solves the equation

\f$-\Delta u+\int u = 1\f$

with homogeneous Dirichlet boundary conditions. This system has a full system matrix which can be written as the sum of of sparse matrix and a rank 1 matrix. The shell matrix concept is used to solve this problem.

The problem is solved in parallel on a non-uniform grid in order to demonstrate all the techniques that are required for this. The grid is fixed, however, i.e. no adaptive mesh refinement is used, so that the example remains simple.

The example is 2d; extension to 3d is straight forward.

C++ include files that we need
        #include <iostream>
        #include <algorithm>
        #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
        #include <cmath>
        
Basic include file needed for the mesh functionality.
        #include "libmesh/libmesh.h"
        #include "libmesh/mesh.h"
        #include "libmesh/mesh_refinement.h"
        #include "libmesh/vtk_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/mesh_generation.h"
        #include "libmesh/sum_shell_matrix.h"
        #include "libmesh/tensor_shell_matrix.h"
        #include "libmesh/sparse_shell_matrix.h"
        #include "libmesh/mesh_refinement.h"
        
        #include "libmesh/getpot.h"
        
This example will solve a linear transient system, so we need to include the \p TransientLinearImplicitSystem definition.
        #include "libmesh/transient_system.h"
        #include "libmesh/linear_implicit_system.h"
        #include "libmesh/vector_value.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 (EquationSystems& es,
        	       const std::string& system_name);
        
Begin the main program. Note that the first statement in the program throws an error if you are in complex number mode, since this example is only intended to work with real numbers.
        int main (int argc, char** argv)
        {
Initialize libMesh.
          LibMeshInit init (argc, argv);
        
        #if !defined(LIBMESH_ENABLE_AMR)
          libmesh_example_assert(false, "--enable-amr");
        #else
          libmesh_example_assert(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
        
Brief message to the user regarding the program name and command line arguments.

          std::cout << "Running: " << argv[0];
        
          for (int i=1; i<argc; i++)
            std::cout << " " << argv[i];
        
          std::cout << std::endl << std::endl;
        
Skip this 2D example if libMesh was compiled as 1D-only.
          libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
        
Create a mesh, with dimension to be overridden later, distributed across the default MPI communicator.
          Mesh mesh(init.comm());
        
Create an equation systems object.
          EquationSystems equation_systems (mesh);
        
          MeshTools::Generation::build_square (mesh,
        				       16,
        				       16,
        				       -1., 1.,
        				       -1., 1.,
        				       QUAD4);
        
          LinearImplicitSystem & system =
            equation_systems.add_system<LinearImplicitSystem>
            ("System");
        
Adds the variable "u" to "System". "u" will be approximated using first-order approximation.
          system.add_variable ("u", FIRST);
        
Also, we need to add two vectors. The tensor matrix v*w^T of these two vectors will be part of the system matrix.
          system.add_vector("v",false);
          system.add_vector("w",false);
        
We need an additional matrix to be used for preconditioning since a shell matrix is not suitable for that.
          system.add_matrix("Preconditioner");
        
Give the system a pointer to the matrix assembly function.
          system.attach_assemble_function (assemble);
        
Initialize the data structures for the equation system.
          equation_systems.init ();
        
Prints information about the system to the screen.
          equation_systems.print_info();
        
          equation_systems.parameters.set<unsigned int>
            ("linear solver maximum iterations") = 250;
          equation_systems.parameters.set<Real>
            ("linear solver tolerance") = TOLERANCE;
        
Refine arbitrarily some elements.
          for(unsigned int i=0; i<2; i++)
            {
              MeshRefinement mesh_refinement(mesh);
              MeshBase::element_iterator       elem_it  = mesh.elements_begin();
              const MeshBase::element_iterator elem_end = mesh.elements_end();
              for (; elem_it != elem_end; ++elem_it)
        	{
        	  Elem* elem = *elem_it;
        	  if(elem->active())
        	    {
        	      if((elem->id()%20)>8)
        		{
        		  elem->set_refinement_flag(Elem::REFINE);
        		}
        	      else
        		{
        		  elem->set_refinement_flag(Elem::DO_NOTHING);
        		}
        	    }
        	  else
        	    {
        	      elem->set_refinement_flag(Elem::INACTIVE);
        	    }
        	}
              mesh_refinement.refine_elements();
              equation_systems.reinit();
            }
        
Prints information about the system to the screen.
          equation_systems.print_info();
        
Before the assemblation of the matrix, we have to clear the two vectors that form the tensor matrix (since this is not performed automatically).
          system.get_vector("v").init(system.n_dofs(), system.n_local_dofs());
          system.get_vector("w").init(system.n_dofs(), system.n_local_dofs());
        
We need a shell matrix to solve. There is currently no way to store the shell matrix in the system. We just create it locally here (a shell matrix does not occupy much memory).
          SumShellMatrix<Number> shellMatrix(system.comm());
          TensorShellMatrix<Number> shellMatrix0(system.get_vector("v"),system.get_vector("w"));
          shellMatrix.matrices.push_back(&shellMatrix0);
          SparseShellMatrix<Number> shellMatrix1(*system.matrix);
          shellMatrix.matrices.push_back(&shellMatrix1);
        
Attach that to the system.
          system.attach_shell_matrix(&shellMatrix);
        
Reset the preconditioning matrix to zero (for the system matrix, the same thing is done automatically).
          system.get_matrix("Preconditioner").zero();
        
Assemble & solve the linear system
          system.solve();
        
Detach the shell matrix from the system since it will go out of scope. Nobody should solve the system outside this function.
          system.detach_shell_matrix();
        
Print a nice message.
          std::cout << "Solved linear system in " << system.n_linear_iterations() << " iterations, residual norm is " << system.final_linear_residual() << "." << std::endl;
        
        #if defined(LIBMESH_HAVE_VTK) && !defined(LIBMESH_ENABLE_PARMESH)
Write result to file.
          VTKIO(mesh).write_equation_systems ("out.pvtu", equation_systems);
        #endif // #ifdef LIBMESH_HAVE_VTK
        
        #endif // #ifndef LIBMESH_ENABLE_AMR
        
          return 0;
        }
        
        
        
This function defines the assembly routine. It is responsible for computing the proper matrix entries for the element stiffness matrices and right-hand sides.
        void assemble (EquationSystems& es,
        	       const std::string& system_name)
        {
        #ifdef LIBMESH_ENABLE_AMR
It is a good idea to make sure we are assembling the proper system.
          libmesh_assert_equal_to (system_name, "System");
        
Get a constant reference to the mesh object.
          const MeshBase& mesh = es.get_mesh();
        
The dimension that we are running
          const unsigned int dim = mesh.mesh_dimension();
        
Get a reference to the Convection-Diffusion system object.
          LinearImplicitSystem & system =
            es.get_system<LinearImplicitSystem> ("System");
        
Get the Finite Element type for the first (and only) variable in the system.
          FEType fe_type = system.variable_type(0);
        
Build a Finite Element object of the specified type. Since the \p FEBase::build() member dynamically creates memory we will store the object as an \p AutoPtr. This can be thought of as a pointer that will clean up after itself.
          AutoPtr<FEBase> fe      (FEBase::build(dim, fe_type));
          AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
        
A Gauss quadrature rule for numerical integration. Let the \p FEType object decide what order rule is appropriate.
          QGauss qrule (dim,   fe_type.default_quadrature_order());
          QGauss qface (dim-1, fe_type.default_quadrature_order());
        
Tell the finite element object to use our quadrature rule.
          fe->attach_quadrature_rule      (&qrule);
          fe_face->attach_quadrature_rule (&qface);
        
Here we define some references to cell-specific data that will be used to assemble the linear system. We will start with the element Jacobian * quadrature weight at each integration point.
          const std::vector<Real>& JxW      = fe->get_JxW();
          const std::vector<Real>& JxW_face = fe_face->get_JxW();
        
The element shape functions evaluated at the quadrature points.
          const std::vector<std::vector<Real> >& phi = fe->get_phi();
          const std::vector<std::vector<Real> >& psi = fe_face->get_phi();
        
The element shape function gradients evaluated at the quadrature points.
          const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
        
The XY locations of the quadrature points used for face integration const std::vector& qface_points = fe_face->get_xyz();

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;
        
Analogous data structures for thw two vectors v and w that form the tensor shell matrix.
          DenseVector<Number> Ve;
          DenseVector<Number> We;
        
This vector will hold the degree of freedom indices for the element. These define where in the global system the element degrees of freedom get mapped.
          std::vector<dof_id_type> dof_indices;
        
Now we will loop over all the elements in the mesh that live on the local processor. We will compute the element matrix and right-hand-side contribution. Since the mesh will be refined we want to 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);
        
Compute the element-specific data for the current element. This involves computing the location of the quadrature points (q_point) and the shape functions (phi, dphi) for the current element.
              fe->reinit (elem);
        
Zero the element matrix and right-hand side before summing them. We use the resize member here because the number of degrees of freedom might have changed from the last element. Note that this will be the case if the element type is different (i.e. the last element was a triangle, now we are on a quadrilateral).
              Ke.resize (dof_indices.size(),
                         dof_indices.size());
        
              Fe.resize (dof_indices.size());
              Ve.resize (dof_indices.size());
              We.resize (dof_indices.size());
        
Now we will build the element matrix and right-hand-side. Constructing the RHS requires the solution and its gradient from the previous timestep. This myst be calculated at each quadrature point by summing the solution degree-of-freedom values by the appropriate weight functions.
              for (unsigned int qp=0; qp<qrule.n_points(); qp++)
                {
Now compute the element matrix and RHS contributions.
                  for (unsigned int i=0; i<phi.size(); i++)
                    {
The RHS contribution
                      Fe(i) += JxW[qp]*(
                                        phi[i][qp]
                                        );
        
                      for (unsigned int j=0; j<phi.size(); j++)
                        {
The matrix contribution
                          Ke(i,j) += JxW[qp]*(
Stiffness matrix
                                              (dphi[i][qp]*dphi[j][qp])
                                              );
                        }
        
V and W are the same for this example.
                      Ve(i) += JxW[qp]*(
                                        phi[i][qp]
                                        );
                      We(i) += JxW[qp]*(
                                        phi[i][qp]
                                        );
                    }
                }
        
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 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.
              {
The penalty value.
                const Real penalty = 1.e10;
        
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)
                    {
                      fe_face->reinit(elem,s);
        
                      for (unsigned int qp=0; qp<qface.n_points(); qp++)
                        {
Matrix contribution
                          for (unsigned int i=0; i<psi.size(); i++)
                            for (unsigned int j=0; j<psi.size(); j++)
                              Ke(i,j) += penalty*JxW_face[qp]*psi[i][qp]*psi[j][qp];
                        }
                    }
              }
        
        
We have now built the element matrix and RHS vector in terms of the element degrees of freedom. However, it is possible that some of the element DOFs are constrained to enforce solution continuity, i.e. they are not really "free". We need to constrain those DOFs in terms of non-constrained DOFs to ensure a continuous solution. The \p DofMap::constrain_element_matrix_and_vector() method does just that.

However, constraining both the sparse matrix (and right hand side) plus the rank 1 matrix is tricky. The dof_indices vector has to be backuped for that because the constraining functions modify it.

              std::vector<dof_id_type> dof_indices_backup(dof_indices);
              dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
              dof_indices = dof_indices_backup;
              dof_map.constrain_element_dyad_matrix(Ve,We,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 SparseMatrix::add_matrix() and \p NumericVector::add_vector() members do this for us.
              system.matrix->add_matrix (Ke, dof_indices);
              system.get_matrix("Preconditioner").add_matrix (Ke, dof_indices);
              system.rhs->add_vector    (Fe, dof_indices);
              system.get_vector("v").add_vector(Ve,dof_indices);
              system.get_vector("w").add_vector(We,dof_indices);
            }
Finished computing the sytem matrix and right-hand side.

Matrices and vectors must be closed manually. This is necessary because the matrix is not directly used as the system matrix (in which case the solver closes it) but as a part of a shell matrix.
          system.matrix->close();
          system.get_matrix("Preconditioner").close();
          system.rhs->close();
          system.get_vector("v").close();
          system.get_vector("w").close();
        
        #endif // #ifdef LIBMESH_ENABLE_AMR
        }



The source file miscellaneous_ex4.C without comments:

 
  
  #include <iostream>
  #include <algorithm>
  #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
  #include <cmath>
  
  #include "libmesh/libmesh.h"
  #include "libmesh/mesh.h"
  #include "libmesh/mesh_refinement.h"
  #include "libmesh/vtk_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/mesh_generation.h"
  #include "libmesh/sum_shell_matrix.h"
  #include "libmesh/tensor_shell_matrix.h"
  #include "libmesh/sparse_shell_matrix.h"
  #include "libmesh/mesh_refinement.h"
  
  #include "libmesh/getpot.h"
  
  #include "libmesh/transient_system.h"
  #include "libmesh/linear_implicit_system.h"
  #include "libmesh/vector_value.h"
  
  #include "libmesh/elem.h"
  
  using namespace libMesh;
  
  void assemble (EquationSystems& es,
  	       const std::string& system_name);
  
  int main (int argc, char** argv)
  {
    LibMeshInit init (argc, argv);
  
  #if !defined(LIBMESH_ENABLE_AMR)
    libmesh_example_assert(false, "--enable-amr");
  #else
    libmesh_example_assert(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
  
  
    std::cout << "Running: " << argv[0];
  
    for (int i=1; i<argc; i++)
      std::cout << " " << argv[i];
  
    std::cout << std::endl << std::endl;
  
    libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
  
    Mesh mesh(init.comm());
  
    EquationSystems equation_systems (mesh);
  
    MeshTools::Generation::build_square (mesh,
  				       16,
  				       16,
  				       -1., 1.,
  				       -1., 1.,
  				       QUAD4);
  
    LinearImplicitSystem & system =
      equation_systems.add_system<LinearImplicitSystem>
      ("System");
  
    system.add_variable ("u", FIRST);
  
    system.add_vector("v",false);
    system.add_vector("w",false);
  
    system.add_matrix("Preconditioner");
  
    system.attach_assemble_function (assemble);
  
    equation_systems.init ();
  
    equation_systems.print_info();
  
    equation_systems.parameters.set<unsigned int>
      ("linear solver maximum iterations") = 250;
    equation_systems.parameters.set<Real>
      ("linear solver tolerance") = TOLERANCE;
  
    for(unsigned int i=0; i<2; i++)
      {
        MeshRefinement mesh_refinement(mesh);
        MeshBase::element_iterator       elem_it  = mesh.elements_begin();
        const MeshBase::element_iterator elem_end = mesh.elements_end();
        for (; elem_it != elem_end; ++elem_it)
  	{
  	  Elem* elem = *elem_it;
  	  if(elem->active())
  	    {
  	      if((elem->id()%20)>8)
  		{
  		  elem->set_refinement_flag(Elem::REFINE);
  		}
  	      else
  		{
  		  elem->set_refinement_flag(Elem::DO_NOTHING);
  		}
  	    }
  	  else
  	    {
  	      elem->set_refinement_flag(Elem::INACTIVE);
  	    }
  	}
        mesh_refinement.refine_elements();
        equation_systems.reinit();
      }
  
    equation_systems.print_info();
  
    system.get_vector("v").init(system.n_dofs(), system.n_local_dofs());
    system.get_vector("w").init(system.n_dofs(), system.n_local_dofs());
  
    SumShellMatrix<Number> shellMatrix(system.comm());
    TensorShellMatrix<Number> shellMatrix0(system.get_vector("v"),system.get_vector("w"));
    shellMatrix.matrices.push_back(&shellMatrix0);
    SparseShellMatrix<Number> shellMatrix1(*system.matrix);
    shellMatrix.matrices.push_back(&shellMatrix1);
  
    system.attach_shell_matrix(&shellMatrix);
  
    system.get_matrix("Preconditioner").zero();
  
    system.solve();
  
    system.detach_shell_matrix();
  
    std::cout << "Solved linear system in " << system.n_linear_iterations() << " iterations, residual norm is " << system.final_linear_residual() << "." << std::endl;
  
  #if defined(LIBMESH_HAVE_VTK) && !defined(LIBMESH_ENABLE_PARMESH)
    VTKIO(mesh).write_equation_systems ("out.pvtu", equation_systems);
  #endif // #ifdef LIBMESH_HAVE_VTK
  
  #endif // #ifndef LIBMESH_ENABLE_AMR
  
    return 0;
  }
  
  
  
  void assemble (EquationSystems& es,
  	       const std::string& system_name)
  {
  #ifdef LIBMESH_ENABLE_AMR
    libmesh_assert_equal_to (system_name, "System");
  
    const MeshBase& mesh = es.get_mesh();
  
    const unsigned int dim = mesh.mesh_dimension();
  
    LinearImplicitSystem & system =
      es.get_system<LinearImplicitSystem> ("System");
  
    FEType fe_type = system.variable_type(0);
  
    AutoPtr<FEBase> fe      (FEBase::build(dim, fe_type));
    AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
  
    QGauss qrule (dim,   fe_type.default_quadrature_order());
    QGauss qface (dim-1, fe_type.default_quadrature_order());
  
    fe->attach_quadrature_rule      (&qrule);
    fe_face->attach_quadrature_rule (&qface);
  
    const std::vector<Real>& JxW      = fe->get_JxW();
    const std::vector<Real>& JxW_face = fe_face->get_JxW();
  
    const std::vector<std::vector<Real> >& phi = fe->get_phi();
    const std::vector<std::vector<Real> >& psi = fe_face->get_phi();
  
    const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
  
  
    const DofMap& dof_map = system.get_dof_map();
  
    DenseMatrix<Number> Ke;
    DenseVector<Number> Fe;
  
    DenseVector<Number> Ve;
    DenseVector<Number> We;
  
    std::vector<dof_id_type> dof_indices;
  
    MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
    const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
  
    for ( ; el != end_el; ++el)
      {
        const Elem* elem = *el;
  
        dof_map.dof_indices (elem, dof_indices);
  
        fe->reinit (elem);
  
        Ke.resize (dof_indices.size(),
                   dof_indices.size());
  
        Fe.resize (dof_indices.size());
        Ve.resize (dof_indices.size());
        We.resize (dof_indices.size());
  
        for (unsigned int qp=0; qp<qrule.n_points(); qp++)
          {
            for (unsigned int i=0; i<phi.size(); i++)
              {
                Fe(i) += JxW[qp]*(
                                  phi[i][qp]
                                  );
  
                for (unsigned int j=0; j<phi.size(); j++)
                  {
                    Ke(i,j) += JxW[qp]*(
  				      (dphi[i][qp]*dphi[j][qp])
                                        );
                  }
  
                Ve(i) += JxW[qp]*(
                                  phi[i][qp]
                                  );
                We(i) += JxW[qp]*(
                                  phi[i][qp]
                                  );
              }
          }
  
        {
          const Real penalty = 1.e10;
  
          for (unsigned int s=0; s<elem->n_sides(); s++)
            if (elem->neighbor(s) == NULL)
              {
                fe_face->reinit(elem,s);
  
                for (unsigned int qp=0; qp<qface.n_points(); qp++)
                  {
                    for (unsigned int i=0; i<psi.size(); i++)
                      for (unsigned int j=0; j<psi.size(); j++)
                        Ke(i,j) += penalty*JxW_face[qp]*psi[i][qp]*psi[j][qp];
                  }
              }
        }
  
  
  
  
        std::vector<dof_id_type> dof_indices_backup(dof_indices);
        dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
        dof_indices = dof_indices_backup;
        dof_map.constrain_element_dyad_matrix(Ve,We,dof_indices);
  
        system.matrix->add_matrix (Ke, dof_indices);
        system.get_matrix("Preconditioner").add_matrix (Ke, dof_indices);
        system.rhs->add_vector    (Fe, dof_indices);
        system.get_vector("v").add_vector(Ve,dof_indices);
        system.get_vector("w").add_vector(We,dof_indices);
      }
  
    system.matrix->close();
    system.get_matrix("Preconditioner").close();
    system.rhs->close();
    system.get_vector("v").close();
    system.get_vector("w").close();
  
  #endif // #ifdef LIBMESH_ENABLE_AMR
  }



The console output of the program:

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

 EquationSystems
  n_systems()=1
   System #0, "System"
    Type "LinearImplicit"
    Variables="u" 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="FIRST", "THIRD" 
    n_dofs()=289
    n_local_dofs()=82
    n_constrained_dofs()=0
    n_local_constrained_dofs()=0
    n_vectors()=3
    n_matrices()=2
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 7.99308
      Average Off-Processor Bandwidth <= 0.719723
      Maximum  On-Processor Bandwidth <= 11
      Maximum Off-Processor Bandwidth <= 7
    DofMap Constraints
      Number of DoF Constraints = 0
      Number of Node Constraints = 0

 EquationSystems
  n_systems()=1
   System #0, "System"
    Type "LinearImplicit"
    Variables="u" 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="FIRST", "THIRD" 
    n_dofs()=2093
    n_local_dofs()=546
    n_constrained_dofs()=458
    n_local_constrained_dofs()=109
    n_vectors()=3
    n_matrices()=2
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 9.41806
      Average Off-Processor Bandwidth <= 0.300048
      Maximum  On-Processor Bandwidth <= 21
      Maximum Off-Processor Bandwidth <= 8
    DofMap Constraints
      Number of DoF Constraints = 458
      Average DoF Constraint Length= 2
      Number of Node Constraints = 909
      Maximum Node Constraint Length= 5
      Average Node Constraint Length= 2.51155

Solved linear system in 26 iterations, residual norm is 1.11367e-06.
*** Warning, This code is untested, experimental, or likely to see future API changes: ../src/mesh/vtk_io.C, line 358, compiled Apr 19 2013 at 11:42:41 ***

 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:51:02 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.342003, Active time=0.256288                                                 |
 ----------------------------------------------------------------------------------------------------------------
| 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()     3         0.0038      0.001262    0.0047      0.001565    1.48     1.83     |
|   build_constraint_matrix()        896       0.0023      0.000003    0.0023      0.000003    0.89     0.89     |
|   build_sparsity()                 3         0.0035      0.001182    0.0081      0.002686    1.38     3.14     |
|   cnstrn_elem_dyad_mat()           448       0.0008      0.000002    0.0008      0.000002    0.33     0.33     |
|   cnstrn_elem_mat_vec()            448       0.0081      0.000018    0.0081      0.000018    3.17     3.17     |
|   create_dof_constraints()         3         0.0129      0.004295    0.0253      0.008441    5.03     9.88     |
|   distribute_dofs()                3         0.0082      0.002721    0.0222      0.007392    3.19     8.65     |
|   dof_indices()                    4672      0.0200      0.000004    0.0200      0.000004    7.79     7.79     |
|   enforce_constraints_exactly()    2         0.0005      0.000237    0.0005      0.000237    0.18     0.18     |
|   old_dof_indices()                1234      0.0062      0.000005    0.0062      0.000005    2.42     2.42     |
|   prepare_send_list()              3         0.0000      0.000006    0.0000      0.000006    0.01     0.01     |
|   reinit()                         3         0.0125      0.004178    0.0125      0.004178    4.89     4.89     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          1         0.0014      0.001442    0.0042      0.004216    0.56     1.65     |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        488       0.0016      0.000003    0.0016      0.000003    0.64     0.64     |
|   init_shape_functions()           41        0.0001      0.000002    0.0001      0.000002    0.03     0.03     |
|   inverse_map()                    3532      0.0065      0.000002    0.0065      0.000002    2.54     2.54     |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             488       0.0014      0.000003    0.0014      0.000003    0.55     0.55     |
|   compute_face_map()               40        0.0003      0.000007    0.0006      0.000014    0.10     0.23     |
|   init_face_shape_functions()      1         0.0000      0.000006    0.0000      0.000006    0.00     0.00     |
|   init_reference_to_physical_map() 41        0.0002      0.000005    0.0002      0.000005    0.08     0.08     |
|                                                                                                                |
| LocationMap                                                                                                    |
|   find()                           6156      0.0078      0.000001    0.0078      0.000001    3.06     3.06     |
|   init()                           4         0.0018      0.000446    0.0018      0.000446    0.70     0.70     |
|                                                                                                                |
| Mesh                                                                                                           |
|   contract()                       2         0.0005      0.000230    0.0008      0.000411    0.18     0.32     |
|   find_neighbors()                 3         0.0109      0.003629    0.0114      0.003801    4.25     4.45     |
|   renumber_nodes_and_elem()        8         0.0012      0.000144    0.0012      0.000144    0.45     0.45     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   compute_hilbert_indices()        4         0.0122      0.003057    0.0122      0.003057    4.77     4.77     |
|   find_global_indices()            4         0.0020      0.000498    0.0159      0.003983    0.78     6.22     |
|   parallel_sort()                  4         0.0006      0.000146    0.0011      0.000284    0.23     0.44     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         1         0.0076      0.007561    0.0119      0.011867    2.95     4.63     |
|                                                                                                                |
| MeshRefinement                                                                                                 |
|   _coarsen_elements()              2         0.0004      0.000212    0.0004      0.000219    0.17     0.17     |
|   _refine_elements()               4         0.0116      0.002905    0.0325      0.008130    4.53     12.69    |
|   add_point()                      6156      0.0116      0.000002    0.0201      0.000003    4.51     7.83     |
|   make_coarsening_compatible()     4         0.0044      0.001109    0.0044      0.001109    1.73     1.73     |
|   make_flags_parallel_consistent() 6         0.0050      0.000833    0.0069      0.001154    1.95     2.70     |
|   make_refinement_compatible()     9         0.0010      0.000107    0.0010      0.000113    0.38     0.40     |
|                                                                                                                |
| MeshTools::Generation                                                                                          |
|   build_cube()                     1         0.0003      0.000274    0.0003      0.000274    0.11     0.11     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      3         0.0287      0.009573    0.0432      0.014393    11.21    16.85    |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      17        0.0005      0.000030    0.0006      0.000033    0.20     0.22     |
|   max(bool)                        19        0.0003      0.000016    0.0003      0.000016    0.12     0.12     |
|   max(scalar)                      465       0.0020      0.000004    0.0020      0.000004    0.78     0.78     |
|   max(vector)                      112       0.0008      0.000007    0.0023      0.000021    0.32     0.91     |
|   min(bool)                        578       0.0023      0.000004    0.0023      0.000004    0.91     0.91     |
|   min(scalar)                      451       0.0038      0.000009    0.0038      0.000009    1.50     1.50     |
|   min(vector)                      112       0.0009      0.000008    0.0024      0.000021    0.36     0.92     |
|   probe()                          168       0.0007      0.000004    0.0007      0.000004    0.29     0.29     |
|   receive()                        168       0.0006      0.000003    0.0013      0.000008    0.23     0.52     |
|   send()                           168       0.0003      0.000002    0.0003      0.000002    0.13     0.13     |
|   send_receive()                   176       0.0008      0.000005    0.0027      0.000015    0.31     1.06     |
|   sum()                            39        0.0008      0.000019    0.0010      0.000025    0.30     0.38     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           168       0.0002      0.000001    0.0002      0.000001    0.07     0.07     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         3         0.0030      0.001002    0.0036      0.001214    1.17     1.42     |
|   set_parent_processor_ids()       3         0.0013      0.000426    0.0013      0.000426    0.50     0.50     |
|                                                                                                                |
| PetscLinearSolver                                                                                              |
|   solve()                          1         0.0185      0.018520    0.0185      0.018520    7.23     7.23     |
|                                                                                                                |
| ProjectVector                                                                                                  |
|   operator()                       2         0.0029      0.001444    0.0112      0.005609    1.13     4.38     |
|                                                                                                                |
| System                                                                                                         |
|   assemble()                       1         0.0162      0.016213    0.0344      0.034434    6.33     13.44    |
|   project_vector()                 2         0.0024      0.001191    0.0175      0.008746    0.93     6.83     |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            27374     0.2563                                          100.00            |
 ----------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example miscellaneous_ex4:
*  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/miscellaneous/miscellaneous_ex4'

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

Hosted By:
SourceForge.net Logo