The source file laplace_exact_solution.h with comments: 

        #include "libmesh/libmesh_common.h"
        
        using namespace libMesh;
        
        #ifndef __laplace_exact_solution_h__
        #define __laplace_exact_solution_h__
        
        class LaplaceExactSolution
        {
        public:
          LaplaceExactSolution(){}
        
          ~LaplaceExactSolution(){}
          
          Real operator()( unsigned int component, 
        		   Real x, Real y, Real z = 0.0)
          {
            const Real hp = 0.5*pi;
        
            switch(component)
            {
            case 0:
              return cos(hp*x)*sin(hp*y)*cos(hp*z);
        
            case 1:
              return sin(hp*x)*cos(hp*y)*cos(hp*z);
        
            case 2:
              return sin(hp*x)*cos(hp*y)*sin(hp*z);
        
            default:
              libmesh_error();
            }
          }
        };
        
        
        class LaplaceExactGradient
        {
        public:
          LaplaceExactGradient(){}
        
          ~LaplaceExactGradient(){}
          
          RealGradient operator()( unsigned int component, 
        			   Real x, Real y, Real z = 0.0)
          {
            const Real hp = 0.5*pi;
        
            switch(component)
            {
            case 0:
              return RealGradient( -hp*sin(hp*x)*sin(hp*y)*cos(hp*z),
        			   cos(hp*x)*(hp)*cos(hp*y)*cos(hp*z),
        			   cos(hp*x)*sin(hp*y)*(-hp)*sin(hp*z) );
        
            case 1:
              return RealGradient( hp*cos(hp*x)*cos(hp*y)*cos(hp*z),
        			   sin(hp*x)*(-hp)*sin(hp*y)*cos(hp*z),
        			   sin(hp*x)*cos(hp*y)*(-hp)*sin(hp*z) );
        
            case 2:
              return RealGradient( hp*cos(hp*x)*cos(hp*y)*sin(hp*z),
        			   sin(hp*x)*(-hp)*sin(hp*y)*sin(hp*z),
        			   sin(hp*x)*cos(hp*y)*(hp)*cos(hp*z) );
        
            default:
              libmesh_error();
            }
          }
        };
        
        #endif // __laplace_exact_solution_h__



The source file laplace_system.h with comments: 

        #include "libmesh/fem_system.h"
        #include "libmesh/vector_value.h"
        #include "libmesh/tensor_value.h"
        #include "libmesh/dirichlet_boundaries.h"
        
        #include "solution_function.h"
        
        using namespace libMesh;
        
        #ifndef __laplace_system_h__
        #define __laplace_system_h__
FEMSystem, TimeSolver and NewtonSolver will handle most tasks, but we must specify element residuals 
        class LaplaceSystem : public FEMSystem
        {
        
        public:
Constructor 
          LaplaceSystem( EquationSystems& es,
        		 const std::string& name,
        		 const unsigned int number);
System initialization 
          virtual void init_data ();
Context initialization 
          virtual void init_context(DiffContext &context);
Element residual and jacobian calculations Time dependent parts 
          virtual bool element_time_derivative (bool request_jacobian,
                                                DiffContext& context);
Constraint part 
          /*
          virtual bool side_constraint (bool request_jacobian,
        				DiffContext& context);
          */
        
        protected:
Indices for each variable; 
          unsigned int u_var;
        
          void init_dirichlet_bcs();
Returns the value of a forcing function at point p. 
          RealGradient forcing(const Point& p);
          
          LaplaceExactSolution exact_solution;
        };
        
        #endif //__laplace_system_h__



The source file solution_function.h with comments: 

        #include "libmesh/function_base.h"
        #include "laplace_exact_solution.h"
        
        using namespace libMesh;
        
        #ifndef __solution_function_h__
        #define __solution_function_h__
        
        class SolutionFunction : public FunctionBase<Number>
        {
        public:
        
          SolutionFunction( const unsigned int u_var )
          : _u_var(u_var) {}
          ~SolutionFunction( ){}
        
          virtual Number operator() (const Point&, const Real = 0)
            { libmesh_not_implemented(); }
        
          virtual void operator() (const Point& p,
                                   const Real,
                                   DenseVector<Number>& output)
          {
            output.zero();
            const Real x=p(0), y=p(1), z=p(2);
libMesh assumes each component of the vector-valued variable is stored contiguously. 
            output(_u_var)   = soln( 0, x, y, z );
            output(_u_var+1) = soln( 1, x, y, z );
            output(_u_var+2) = soln( 2, x, y, z );
          }
        
          virtual Number component( unsigned int component_in, const Point& p,
        			    const Real )
          {
            const Real x=p(0), y=p(1), z=p(2);
            return soln( component_in, x, y, z );
          }
        
          virtual AutoPtr<FunctionBase<Number> > clone() const
          { return AutoPtr<FunctionBase<Number> > (new SolutionFunction(_u_var)); }
        
        private:
        
          const unsigned int _u_var;
          LaplaceExactSolution soln;
        };
FIXME: PB: We ought to be able to merge the above class with this one through templating, but I'm being lazy. 
        class SolutionGradient : public FunctionBase<Gradient>
        {
        public:
        
          SolutionGradient( const unsigned int u_var )
          : _u_var(u_var) {}
          ~SolutionGradient( ){}
        
          virtual Gradient operator() (const Point&, const Real = 0)
            { libmesh_not_implemented(); }
        
          virtual void operator() (const Point& p,
                                   const Real,
                                   DenseVector<Gradient>& output)
          {
            output.zero();
            const Real x=p(0), y=p(1), z=p(2);
            output(_u_var)   = soln( 0, x, y, z );
            output(_u_var+1) = soln( 1, x, y, z );
            output(_u_var+2) = soln( 2, x, y, z );
          }
        
          virtual Gradient component( unsigned int component_in, const Point& p,
        			    const Real )
          {
            const Real x=p(0), y=p(1), z=p(2);
            return soln( component_in, x, y, z );
          }
        
          virtual AutoPtr<FunctionBase<Gradient> > clone() const
          { return AutoPtr<FunctionBase<Gradient> > (new SolutionGradient(_u_var)); }
        
        private:
        
          const unsigned int _u_var;
          LaplaceExactGradient soln;
        };
        
        #endif // __solution_function_h__



The source file laplace_system.C with comments: 

        #include "libmesh/getpot.h"
        
        #include "laplace_system.h"
        
        #include "libmesh/boundary_info.h"
        #include "libmesh/dof_map.h"
        #include "libmesh/fe_base.h"
        #include "libmesh/fe_interface.h"
        #include "libmesh/fem_context.h"
        #include "libmesh/mesh.h"
        #include "libmesh/quadrature.h"
        #include "libmesh/string_to_enum.h"
Bring in everything from the libMesh namespace 
        using namespace libMesh;
        
        LaplaceSystem::LaplaceSystem( EquationSystems& es,
        			      const std::string& name_in,
        			      const unsigned int number_in)
          : FEMSystem(es, name_in, number_in)
        {
          return;
        }
        
        void LaplaceSystem::init_data ()
        {
Check the input file for Reynolds number, application type, approximation type 
          GetPot infile("laplace.in");
Add the solution variable 
          u_var = this->add_variable ("u", FIRST, LAGRANGE_VEC);
        
          this->time_evolving(u_var);
Useful debugging options Set verify_analytic_jacobians to 1e-6 to use 
          this->verify_analytic_jacobians = infile("verify_analytic_jacobians", 0.);
          this->print_jacobians = infile("print_jacobians", false);
          this->print_element_jacobians = infile("print_element_jacobians", false);
        
          this->init_dirichlet_bcs();
Do the parent's initialization after variables and boundary constraints are defined 
          FEMSystem::init_data();
        }
        
        void LaplaceSystem::init_dirichlet_bcs()
        {
          const boundary_id_type all_ids[6] = {0,1,2,3,4,5};
          std::set<boundary_id_type> boundary_ids( all_ids, all_ids+6 );
        
          std::vector<unsigned int> vars;
          vars.push_back( u_var );
Note that for vector-valued variables, it is assumed each component is stored contiguously. For 2-D elements in 3-D space, only two components should be returned. 
          SolutionFunction func( u_var );
        
          this->get_dof_map().add_dirichlet_boundary( libMesh::DirichletBoundary( boundary_ids, vars, &func ) );
          return;
        }
        
        void LaplaceSystem::init_context(DiffContext &context)
        {
          FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
Get finite element object 
          FEGenericBase<RealGradient>* fe;
          c.get_element_fe<RealGradient>( u_var, fe );
We should prerequest all the data we will need to build the linear system. 
          fe->get_JxW();
          fe->get_phi();
          fe->get_dphi();
          fe->get_xyz();
        
          FEGenericBase<RealGradient>* side_fe;
          c.get_side_fe<RealGradient>( u_var, side_fe );
        
          side_fe->get_JxW();
          side_fe->get_phi();
        }
        
        
        bool LaplaceSystem::element_time_derivative (bool request_jacobian,
                                                    DiffContext &context)
        {
          FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
Get finite element object 
          FEGenericBase<RealGradient>* fe = NULL;
          c.get_element_fe<RealGradient>( u_var, fe );
First we get some references to cell-specific data that will be used to assemble the linear system.

Element Jacobian * quadrature weights for interior integration 
          const std::vector<Real> &JxW = fe->get_JxW();
The velocity shape functions at interior quadrature points. 
          const std::vector<std::vector<RealGradient> >& phi = fe->get_phi();
The velocity shape function gradients at interior quadrature points. 
          const std::vector<std::vector<RealTensor> >& grad_phi = fe->get_dphi();
         
          const std::vector<Point>& qpoint = fe->get_xyz();
The number of local degrees of freedom in each variable 
          const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
        
          DenseSubMatrix<Number> &Kuu = *c.elem_subjacobians[u_var][u_var];
        
          DenseSubVector<Number> &Fu = *c.elem_subresiduals[u_var];
Now we will build the element Jacobian and residual. Constructing the residual requires the solution and its gradient from the previous timestep. This must be calculated at each quadrature point by summing the solution degree-of-freedom values by the appropriate weight functions. 
          const unsigned int n_qpoints = (c.get_element_qrule())->n_points();
        
          for (unsigned int qp=0; qp != n_qpoints; qp++)
            {
              Tensor grad_u;
              
              c.interior_gradient( u_var, qp, grad_u );
Value of the forcing function at this quadrature point 
              RealGradient f = this->forcing(qpoint[qp]);
First, an i-loop over the velocity degrees of freedom. We know that n_u_dofs == n_v_dofs so we can compute contributions for both at the same time. 
              for (unsigned int i=0; i != n_u_dofs; i++)
                {
                  Fu(i) += ( grad_u.contract(grad_phi[i][qp]) - f*phi[i][qp] )*JxW[qp];
        
                  if (request_jacobian)
                    {
Matrix contributions for the uu and vv couplings. 
                      for (unsigned int j=0; j != n_u_dofs; j++)
                        {
                          Kuu(i,j) += grad_phi[j][qp].contract(grad_phi[i][qp])*JxW[qp];
                         
        		}
        	    }
                }
            } // end of the quadrature point qp-loop
          
          return request_jacobian;
        }
        
        /*
        bool LaplaceSystem::side_constraint (bool request_jacobian,
        				     DiffContext &context)
        {
          FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
Get finite element object 
          FEGenericBase<RealGradient>* side_fe = NULL;
          c.get_side_fe<RealGradient>( u_var, side_fe );
First we get some references to cell-specific data that will be used to assemble the linear system.

Element Jacobian * quadrature weights for interior integration 
          const std::vector<Real> &JxW = side_fe->get_JxW();
The velocity shape functions at interior quadrature points. 
          const std::vector<std::vector<RealGradient> >& phi = side_fe->get_phi();
The number of local degrees of freedom in each variable 
          const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
        
          const std::vector<Point>& qpoint = side_fe->get_xyz();
The penalty value. \frac{1}{\epsilon} in the discussion above. 
          const Real penalty = 1.e10;
        
          DenseSubMatrix<Number> &Kuu = *c.elem_subjacobians[u_var][u_var];
          DenseSubVector<Number> &Fu = *c.elem_subresiduals[u_var];
        
          const unsigned int n_qpoints = (c.get_side_qrule())->n_points();
        
          for (unsigned int qp=0; qp != n_qpoints; qp++)
            {
              Gradient u;
              c.side_value( u_var, qp, u );
         
              Gradient u_exact( this->exact_solution( 0, qpoint[qp](0), qpoint[qp](1) ),
        			this->exact_solution( 1, qpoint[qp](0), qpoint[qp](1) ));
        
              for (unsigned int i=0; i != n_u_dofs; i++)
        	{
        	  Fu(i) += penalty*(u - u_exact)*phi[i][qp]*JxW[qp];
        	  
        	  if (request_jacobian)
        	    {
        	      for (unsigned int j=0; j != n_u_dofs; j++)
        		Kuu(i,j) += penalty*phi[j][qp]*phi[i][qp]*JxW[qp];
        	    }
        	}
            }
        
          return request_jacobian;
        }
        */
        
        RealGradient LaplaceSystem::forcing( const Point& p )
        {
          Real x = p(0); Real y = p(1); Real z = p(2);
          
          const Real eps = 1.e-3;
        
          const Real fx = -(exact_solution(0,x,y,z-eps) +
        		    exact_solution(0,x,y,z+eps) +
        		    exact_solution(0,x,y-eps,z) +
        		    exact_solution(0,x,y+eps,z) +
        		    exact_solution(0,x-eps,y,z) +
        		    exact_solution(0,x+eps,y,z) -
        		    6.*exact_solution(0,x,y,z))/eps/eps;
          
          const Real fy = -(exact_solution(1,x,y,z-eps) +
        		    exact_solution(1,x,y,z+eps) +
        		    exact_solution(1,x,y-eps,z) +
        		    exact_solution(1,x,y+eps,z) +
        		    exact_solution(1,x-eps,y,z) +
        		    exact_solution(1,x+eps,y,z) -
        		    6.*exact_solution(1,x,y,z))/eps/eps;
        
          const Real fz = -(exact_solution(2,x,y,z-eps) +
        		    exact_solution(2,x,y,z+eps) +
        		    exact_solution(2,x,y-eps,z) +
        		    exact_solution(2,x,y+eps,z) +
        		    exact_solution(2,x-eps,y,z) +
        		    exact_solution(2,x+eps,y,z) -
        		    6.*exact_solution(2,x,y,z))/eps/eps;
          
          return RealGradient( fx, fy, fz );
        }



The source file vector_fe_ex2.C with comments:

FEMSystem Example 1 - Unsteady Navier-Stokes Equations with FEMSystem



This example shows how the transient nonlinear problem from example 13 can be solved using the DifferentiableSystem class framework

Basic include files 
        #include "libmesh/equation_systems.h"
        #include "libmesh/getpot.h"
        #include "libmesh/exodusII_io.h"
        #include "libmesh/mesh.h"
        #include "libmesh/mesh_generation.h"
        #include "libmesh/exact_solution.h"
        #include "libmesh/ucd_io.h"
The systems and solvers we may use 
        #include "laplace_system.h"
        #include "libmesh/diff_solver.h"
        #include "libmesh/steady_solver.h"
Bring in everything from the libMesh namespace 
        using namespace libMesh;
The main program. 
        int main (int argc, char** argv)
        {
Initialize libMesh. 
          LibMeshInit init (argc, argv);
Parse the input file 
          GetPot infile("vector_fe_ex2.in");
Read in parameters from the input file 
          const unsigned int grid_size = infile( "grid_size", 2 );
Skip higher-dimensional examples on a lower-dimensional libMesh build 
          libmesh_example_assert(3 <= LIBMESH_DIM, "2D/3D support");
Create a mesh. 
          Mesh mesh;
Use the MeshTools::Generation mesh generator to create a uniform grid on the square [-1,1]^D. We instruct the mesh generator to build a mesh of 8x8 \p Quad9 elements in 2D, or \p Hex27 elements in 3D. Building these higher-order elements allows us to use higher-order approximation, as in example 3. 
          MeshTools::Generation::build_cube (mesh,
        				     grid_size,
        				     grid_size,
        				     grid_size,
        				     -1., 1.,
        				     -1., 1.,
        				     -1., 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 "Navier-Stokes" and its variables. 
          LaplaceSystem & system = 
            equation_systems.add_system<LaplaceSystem> ("Laplace");
This example only implements the steady-state problem 
          system.time_solver =
            AutoPtr<TimeSolver>(new SteadySolver(system));
Initialize the system 
          equation_systems.init();
And the nonlinear solver options 
          DiffSolver &solver = *(system.time_solver->diff_solver().get());
          solver.quiet = infile("solver_quiet", true);
          solver.verbose = !solver.quiet;
          solver.max_nonlinear_iterations =
            infile("max_nonlinear_iterations", 15);
          solver.relative_step_tolerance =
            infile("relative_step_tolerance", 1.e-3);
          solver.relative_residual_tolerance =
            infile("relative_residual_tolerance", 0.0);
          solver.absolute_residual_tolerance =
            infile("absolute_residual_tolerance", 0.0);
And the linear solver options 
          solver.max_linear_iterations =
            infile("max_linear_iterations", 50000);
          solver.initial_linear_tolerance =
            infile("initial_linear_tolerance", 1.e-3);
Print information about the system to the screen. 
          equation_systems.print_info();
        
          system.solve();
        
          ExactSolution exact_sol( equation_systems );
          
          std::vector<FunctionBase<Number>* > sols;
          std::vector<FunctionBase<Gradient>* > grads;
        
          sols.push_back( new SolutionFunction(system.variable_number("u")) );
          grads.push_back( new SolutionGradient(system.variable_number("u")) );
        
          exact_sol.attach_exact_values(sols);
          exact_sol.attach_exact_derivs(grads);
Use higher quadrature order for more accurate error results 
          int extra_error_quadrature = infile("extra_error_quadrature",2);
          exact_sol.extra_quadrature_order(extra_error_quadrature);
Compute the error. 
          exact_sol.compute_error("Laplace", "u");
Print out the error values 
          std::cout << "L2-Error is: "
        	    << exact_sol.l2_error("Laplace", "u")
        	    << std::endl;
          std::cout << "H1-Error is: "
        	    << exact_sol.h1_error("Laplace", "u")
        	    << std::endl;
          
        #ifdef LIBMESH_HAVE_EXODUS_API
We write the file in the ExodusII format. 
          ExodusII_IO(mesh).write_equation_systems("out.e", equation_systems);
          
        #endif // #ifdef LIBMESH_HAVE_EXODUS_API
        
          UCDIO(mesh).write_equation_systems("out.inp", equation_systems);
All done. 
          return 0;
        }



The source file laplace_exact_solution.h without comments: 

 
  #include "libmesh/libmesh_common.h"
  
  using namespace libMesh;
  
  #ifndef __laplace_exact_solution_h__
  #define __laplace_exact_solution_h__
  
  class LaplaceExactSolution
  {
  public:
    LaplaceExactSolution(){}
  
    ~LaplaceExactSolution(){}
    
    Real operator()( unsigned int component, 
  		   Real x, Real y, Real z = 0.0)
    {
      const Real hp = 0.5*pi;
  
      switch(component)
      {
      case 0:
        return cos(hp*x)*sin(hp*y)*cos(hp*z);
  
      case 1:
        return sin(hp*x)*cos(hp*y)*cos(hp*z);
  
      case 2:
        return sin(hp*x)*cos(hp*y)*sin(hp*z);
  
      default:
        libmesh_error();
      }
    }
  };
  
  
  class LaplaceExactGradient
  {
  public:
    LaplaceExactGradient(){}
  
    ~LaplaceExactGradient(){}
    
    RealGradient operator()( unsigned int component, 
  			   Real x, Real y, Real z = 0.0)
    {
      const Real hp = 0.5*pi;
  
      switch(component)
      {
      case 0:
        return RealGradient( -hp*sin(hp*x)*sin(hp*y)*cos(hp*z),
  			   cos(hp*x)*(hp)*cos(hp*y)*cos(hp*z),
  			   cos(hp*x)*sin(hp*y)*(-hp)*sin(hp*z) );
  
      case 1:
        return RealGradient( hp*cos(hp*x)*cos(hp*y)*cos(hp*z),
  			   sin(hp*x)*(-hp)*sin(hp*y)*cos(hp*z),
  			   sin(hp*x)*cos(hp*y)*(-hp)*sin(hp*z) );
  
      case 2:
        return RealGradient( hp*cos(hp*x)*cos(hp*y)*sin(hp*z),
  			   sin(hp*x)*(-hp)*sin(hp*y)*sin(hp*z),
  			   sin(hp*x)*cos(hp*y)*(hp)*cos(hp*z) );
  
      default:
        libmesh_error();
      }
    }
  };
  
  #endif // __laplace_exact_solution_h__



The source file laplace_system.h without comments: 

 
  #include "libmesh/fem_system.h"
  #include "libmesh/vector_value.h"
  #include "libmesh/tensor_value.h"
  #include "libmesh/dirichlet_boundaries.h"
  
  #include "solution_function.h"
  
  using namespace libMesh;
  
  #ifndef __laplace_system_h__
  #define __laplace_system_h__
  
  class LaplaceSystem : public FEMSystem
  {
  
  public:
  
    LaplaceSystem( EquationSystems& es,
  		 const std::string& name,
  		 const unsigned int number);
    
    virtual void init_data ();
    
    virtual void init_context(DiffContext &context);
    
    virtual bool element_time_derivative (bool request_jacobian,
                                          DiffContext& context);
    
    /*
    virtual bool side_constraint (bool request_jacobian,
  				DiffContext& context);
    */
  
  protected:
    unsigned int u_var;
  
    void init_dirichlet_bcs();
    
    RealGradient forcing(const Point& p);
    
    LaplaceExactSolution exact_solution;
  };
  
  #endif //__laplace_system_h__



The source file solution_function.h without comments: 

 
  #include "libmesh/function_base.h"
  #include "laplace_exact_solution.h"
  
  using namespace libMesh;
  
  #ifndef __solution_function_h__
  #define __solution_function_h__
  
  class SolutionFunction : public FunctionBase<Number>
  {
  public:
  
    SolutionFunction( const unsigned int u_var )
    : _u_var(u_var) {}
    ~SolutionFunction( ){}
  
    virtual Number operator() (const Point&, const Real = 0)
      { libmesh_not_implemented(); }
  
    virtual void operator() (const Point& p,
                             const Real,
                             DenseVector<Number>& output)
    {
      output.zero();
      const Real x=p(0), y=p(1), z=p(2);
      output(_u_var)   = soln( 0, x, y, z );
      output(_u_var+1) = soln( 1, x, y, z );
      output(_u_var+2) = soln( 2, x, y, z );
    }
  
    virtual Number component( unsigned int component_in, const Point& p,
  			    const Real )
    {
      const Real x=p(0), y=p(1), z=p(2);
      return soln( component_in, x, y, z );
    }
  
    virtual AutoPtr<FunctionBase<Number> > clone() const
    { return AutoPtr<FunctionBase<Number> > (new SolutionFunction(_u_var)); }
  
  private:
  
    const unsigned int _u_var;
    LaplaceExactSolution soln;
  };
  
  class SolutionGradient : public FunctionBase<Gradient>
  {
  public:
  
    SolutionGradient( const unsigned int u_var )
    : _u_var(u_var) {}
    ~SolutionGradient( ){}
  
    virtual Gradient operator() (const Point&, const Real = 0)
      { libmesh_not_implemented(); }
  
    virtual void operator() (const Point& p,
                             const Real,
                             DenseVector<Gradient>& output)
    {
      output.zero();
      const Real x=p(0), y=p(1), z=p(2);
      output(_u_var)   = soln( 0, x, y, z );
      output(_u_var+1) = soln( 1, x, y, z );
      output(_u_var+2) = soln( 2, x, y, z );
    }
  
    virtual Gradient component( unsigned int component_in, const Point& p,
  			    const Real )
    {
      const Real x=p(0), y=p(1), z=p(2);
      return soln( component_in, x, y, z );
    }
  
    virtual AutoPtr<FunctionBase<Gradient> > clone() const
    { return AutoPtr<FunctionBase<Gradient> > (new SolutionGradient(_u_var)); }
  
  private:
  
    const unsigned int _u_var;
    LaplaceExactGradient soln;
  };
  
  #endif // __solution_function_h__



The source file laplace_system.C without comments: 

 
  #include "libmesh/getpot.h"
  
  #include "laplace_system.h"
  
  #include "libmesh/boundary_info.h"
  #include "libmesh/dof_map.h"
  #include "libmesh/fe_base.h"
  #include "libmesh/fe_interface.h"
  #include "libmesh/fem_context.h"
  #include "libmesh/mesh.h"
  #include "libmesh/quadrature.h"
  #include "libmesh/string_to_enum.h"
  
  
  using namespace libMesh;
  
  LaplaceSystem::LaplaceSystem( EquationSystems& es,
  			      const std::string& name_in,
  			      const unsigned int number_in)
    : FEMSystem(es, name_in, number_in)
  {
    return;
  }
  
  void LaplaceSystem::init_data ()
  {
    GetPot infile("laplace.in");
  
    u_var = this->add_variable ("u", FIRST, LAGRANGE_VEC);
  
    this->time_evolving(u_var);
  
    this->verify_analytic_jacobians = infile("verify_analytic_jacobians", 0.);
    this->print_jacobians = infile("print_jacobians", false);
    this->print_element_jacobians = infile("print_element_jacobians", false);
  
    this->init_dirichlet_bcs();
  
    FEMSystem::init_data();
  }
  
  void LaplaceSystem::init_dirichlet_bcs()
  {
    const boundary_id_type all_ids[6] = {0,1,2,3,4,5};
    std::set<boundary_id_type> boundary_ids( all_ids, all_ids+6 );
  
    std::vector<unsigned int> vars;
    vars.push_back( u_var );
  
    SolutionFunction func( u_var );
  
    this->get_dof_map().add_dirichlet_boundary( libMesh::DirichletBoundary( boundary_ids, vars, &func ) );
    return;
  }
  
  void LaplaceSystem::init_context(DiffContext &context)
  {
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
  
    FEGenericBase<RealGradient>* fe;
    c.get_element_fe<RealGradient>( u_var, fe );
    
    fe->get_JxW();
    fe->get_phi();
    fe->get_dphi();
    fe->get_xyz();
  
    FEGenericBase<RealGradient>* side_fe;
    c.get_side_fe<RealGradient>( u_var, side_fe );
  
    side_fe->get_JxW();
    side_fe->get_phi();
  }
  
  
  bool LaplaceSystem::element_time_derivative (bool request_jacobian,
                                              DiffContext &context)
  {
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
  
    FEGenericBase<RealGradient>* fe = NULL;
    c.get_element_fe<RealGradient>( u_var, fe );
    
  
    const std::vector<Real> &JxW = fe->get_JxW();
  
    const std::vector<std::vector<RealGradient> >& phi = fe->get_phi();
  
    const std::vector<std::vector<RealTensor> >& grad_phi = fe->get_dphi();
   
    const std::vector<Point>& qpoint = fe->get_xyz();
  
    const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
  
    DenseSubMatrix<Number> &Kuu = *c.elem_subjacobians[u_var][u_var];
  
    DenseSubVector<Number> &Fu = *c.elem_subresiduals[u_var];
        
    const unsigned int n_qpoints = (c.get_element_qrule())->n_points();
  
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        Tensor grad_u;
        
        c.interior_gradient( u_var, qp, grad_u );
  
  
        RealGradient f = this->forcing(qpoint[qp]);
  
        for (unsigned int i=0; i != n_u_dofs; i++)
          {
            Fu(i) += ( grad_u.contract(grad_phi[i][qp]) - f*phi[i][qp] )*JxW[qp];
  
            if (request_jacobian)
              {
                for (unsigned int j=0; j != n_u_dofs; j++)
                  {
                    Kuu(i,j) += grad_phi[j][qp].contract(grad_phi[i][qp])*JxW[qp];
                   
  		}
  	    }
          }
      } // end of the quadrature point qp-loop
    
    return request_jacobian;
  }
  
  /*
  bool LaplaceSystem::side_constraint (bool request_jacobian,
  				     DiffContext &context)
  {
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
  
    FEGenericBase<RealGradient>* side_fe = NULL;
    c.get_side_fe<RealGradient>( u_var, side_fe );
    
  
    const std::vector<Real> &JxW = side_fe->get_JxW();
  
    const std::vector<std::vector<RealGradient> >& phi = side_fe->get_phi();
  
    const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
  
    const std::vector<Point>& qpoint = side_fe->get_xyz();
  
    const Real penalty = 1.e10;
  
    DenseSubMatrix<Number> &Kuu = *c.elem_subjacobians[u_var][u_var];
    DenseSubVector<Number> &Fu = *c.elem_subresiduals[u_var];
  
    const unsigned int n_qpoints = (c.get_side_qrule())->n_points();
  
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        Gradient u;
        c.side_value( u_var, qp, u );
   
        Gradient u_exact( this->exact_solution( 0, qpoint[qp](0), qpoint[qp](1) ),
  			this->exact_solution( 1, qpoint[qp](0), qpoint[qp](1) ));
  
        for (unsigned int i=0; i != n_u_dofs; i++)
  	{
  	  Fu(i) += penalty*(u - u_exact)*phi[i][qp]*JxW[qp];
  	  
  	  if (request_jacobian)
  	    {
  	      for (unsigned int j=0; j != n_u_dofs; j++)
  		Kuu(i,j) += penalty*phi[j][qp]*phi[i][qp]*JxW[qp];
  	    }
  	}
      }
  
    return request_jacobian;
  }
  */
  
  RealGradient LaplaceSystem::forcing( const Point& p )
  {
    Real x = p(0); Real y = p(1); Real z = p(2);
    
    const Real eps = 1.e-3;
  
    const Real fx = -(exact_solution(0,x,y,z-eps) +
  		    exact_solution(0,x,y,z+eps) +
  		    exact_solution(0,x,y-eps,z) +
  		    exact_solution(0,x,y+eps,z) +
  		    exact_solution(0,x-eps,y,z) +
  		    exact_solution(0,x+eps,y,z) -
  		    6.*exact_solution(0,x,y,z))/eps/eps;
    
    const Real fy = -(exact_solution(1,x,y,z-eps) +
  		    exact_solution(1,x,y,z+eps) +
  		    exact_solution(1,x,y-eps,z) +
  		    exact_solution(1,x,y+eps,z) +
  		    exact_solution(1,x-eps,y,z) +
  		    exact_solution(1,x+eps,y,z) -
  		    6.*exact_solution(1,x,y,z))/eps/eps;
  
    const Real fz = -(exact_solution(2,x,y,z-eps) +
  		    exact_solution(2,x,y,z+eps) +
  		    exact_solution(2,x,y-eps,z) +
  		    exact_solution(2,x,y+eps,z) +
  		    exact_solution(2,x-eps,y,z) +
  		    exact_solution(2,x+eps,y,z) -
  		    6.*exact_solution(2,x,y,z))/eps/eps;
    
    return RealGradient( fx, fy, fz );
  }



The source file vector_fe_ex2.C without comments: 

 
  
  #include "libmesh/equation_systems.h"
  #include "libmesh/getpot.h"
  #include "libmesh/exodusII_io.h"
  #include "libmesh/mesh.h"
  #include "libmesh/mesh_generation.h"
  #include "libmesh/exact_solution.h"
  #include "libmesh/ucd_io.h"
  
  #include "laplace_system.h"
  #include "libmesh/diff_solver.h"
  #include "libmesh/steady_solver.h"
  
  using namespace libMesh;
  
  int main (int argc, char** argv)
  {
    LibMeshInit init (argc, argv);
  
    GetPot infile("vector_fe_ex2.in");
  
    const unsigned int grid_size = infile( "grid_size", 2 );
  
    libmesh_example_assert(3 <= LIBMESH_DIM, "2D/3D support");
  
    Mesh mesh;
    
    MeshTools::Generation::build_cube (mesh,
  				     grid_size,
  				     grid_size,
  				     grid_size,
  				     -1., 1.,
  				     -1., 1.,
  				     -1., 1.,
  				     HEX8);
  
    mesh.print_info();
  
    EquationSystems equation_systems (mesh);
  
    LaplaceSystem & system = 
      equation_systems.add_system<LaplaceSystem> ("Laplace");
  
    system.time_solver =
      AutoPtr<TimeSolver>(new SteadySolver(system));
  
    equation_systems.init();
  
    DiffSolver &solver = *(system.time_solver->diff_solver().get());
    solver.quiet = infile("solver_quiet", true);
    solver.verbose = !solver.quiet;
    solver.max_nonlinear_iterations =
      infile("max_nonlinear_iterations", 15);
    solver.relative_step_tolerance =
      infile("relative_step_tolerance", 1.e-3);
    solver.relative_residual_tolerance =
      infile("relative_residual_tolerance", 0.0);
    solver.absolute_residual_tolerance =
      infile("absolute_residual_tolerance", 0.0);
      
    solver.max_linear_iterations =
      infile("max_linear_iterations", 50000);
    solver.initial_linear_tolerance =
      infile("initial_linear_tolerance", 1.e-3);
  
    equation_systems.print_info();
  
    system.solve();
  
    ExactSolution exact_sol( equation_systems );
    
    std::vector<FunctionBase<Number>* > sols;
    std::vector<FunctionBase<Gradient>* > grads;
  
    sols.push_back( new SolutionFunction(system.variable_number("u")) );
    grads.push_back( new SolutionGradient(system.variable_number("u")) );
  
    exact_sol.attach_exact_values(sols);
    exact_sol.attach_exact_derivs(grads);
    
    int extra_error_quadrature = infile("extra_error_quadrature",2);
    exact_sol.extra_quadrature_order(extra_error_quadrature);
  
    exact_sol.compute_error("Laplace", "u");
    
    std::cout << "L2-Error is: "
  	    << exact_sol.l2_error("Laplace", "u")
  	    << std::endl;
    std::cout << "H1-Error is: "
  	    << exact_sol.h1_error("Laplace", "u")
  	    << std::endl;
    
  #ifdef LIBMESH_HAVE_EXODUS_API
    
    ExodusII_IO(mesh).write_equation_systems("out.e", equation_systems);
    
  #endif // #ifdef LIBMESH_HAVE_EXODUS_API
  
    UCDIO(mesh).write_equation_systems("out.inp", equation_systems);
  
    return 0;
  }



The console output of the program: 

***************************************************************
* Running Example vector_fe_ex2:
*  mpirun -np 2 example-devel  -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************
 
 Mesh Information:
  mesh_dimension()=3
  spatial_dimension()=3
  n_nodes()=216
    n_local_nodes()=132
  n_elem()=125
    n_local_elem()=62
    n_active_elem()=125
  n_subdomains()=1
  n_partitions()=2
  n_processors()=2
  n_threads()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System #0, "Laplace"
    Type "Implicit"
    Variables="u" 
    Finite Element Types="LAGRANGE_VEC", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="FIRST", "THIRD" 
    n_dofs()=648
    n_local_dofs()=396
    n_constrained_dofs()=456
    n_local_constrained_dofs()=270
    n_vectors()=1
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 53.6389
      Average Off-Processor Bandwidth <= 7.88889
      Maximum  On-Processor Bandwidth <= 99
      Maximum Off-Processor Bandwidth <= 42
    DofMap Constraints
      Number of DoF Constraints = 456
      Number of Heterogenous Constraints= 456
      Average DoF Constraint Length= 0
      Number of Node Constraints = 0

Assembling the System
Nonlinear Residual: 4.00845
Linear solve starting, tolerance 1e-12
Linear solve finished, step 13, residual 2.10554e-12
Trying full Newton step
  Current Residual: 2.20626e-12
  Nonlinear solver converged, step 0, residual reduction 5.50403e-13 < 1e-12
  Nonlinear solver relative step size 0.723571 > 1e-06
L2-Error is: 0.142039
H1-Error is: 0.902143
************************************************************************************************************************
***             WIDEN YOUR WINDOW TO 120 CHARACTERS.  Use 'enscript -r -fCourier9' to print this document            ***
************************************************************************************************************************

---------------------------------------------- PETSc Performance Summary: ----------------------------------------------

/workspace/libmesh/examples/vector_fe/vector_fe_ex2/.libs/lt-example-devel on a intel-12. named hbar.ices.utexas.edu with 2 processors, by benkirk Tue Feb  5 13:54:13 2013
Using Petsc Release Version 3.3.0, Patch 2, Fri Jul 13 15:42:00 CDT 2012 

                         Max       Max/Min        Avg      Total 
Time (sec):           9.525e-01      1.00000   9.525e-01
Objects:              6.700e+01      1.00000   6.700e+01
Flops:                3.276e+06      2.40345   2.320e+06  4.639e+06
Flops/sec:            3.439e+06      2.40345   2.435e+06  4.870e+06
MPI Messages:         3.450e+01      1.00000   3.450e+01  6.900e+01
MPI Message Lengths:  1.175e+05      1.00000   3.406e+03  2.350e+05
MPI Reductions:       1.330e+02      1.00000

Flop counting convention: 1 flop = 1 real number operation of type (multiply/divide/add/subtract)
                            e.g., VecAXPY() for real vectors of length N --> 2N flops
                            and VecAXPY() for complex vectors of length N --> 8N flops

Summary of Stages:   ----- Time ------  ----- Flops -----  --- Messages ---  -- Message Lengths --  -- Reductions --
                        Avg     %Total     Avg     %Total   counts   %Total     Avg         %Total   counts   %Total 
 0:      Main Stage: 9.5250e-01 100.0%  4.6391e+06 100.0%  6.900e+01 100.0%  3.406e+03      100.0%  1.320e+02  99.2% 

------------------------------------------------------------------------------------------------------------------------
See the 'Profiling' chapter of the users' manual for details on interpreting output.
Phase summary info:
   Count: number of times phase was executed
   Time and Flops: Max - maximum over all processors
                   Ratio - ratio of maximum to minimum over all processors
   Mess: number of messages sent
   Avg. len: average message length
   Reduct: number of global reductions
   Global: entire computation
   Stage: stages of a computation. Set stages with PetscLogStagePush() and PetscLogStagePop().
      %T - percent time in this phase         %f - percent flops in this phase
      %M - percent messages in this phase     %L - percent message lengths in this phase
      %R - percent reductions in this phase
   Total Mflop/s: 10e-6 * (sum of flops over all processors)/(max time over all processors)
------------------------------------------------------------------------------------------------------------------------
Event                Count      Time (sec)     Flops                             --- Global ---  --- Stage ---   Total
                   Max Ratio  Max     Ratio   Max  Ratio  Mess   Avg len Reduct  %T %f %M %L %R  %T %f %M %L %R Mflop/s
------------------------------------------------------------------------------------------------------------------------

--- Event Stage 0: Main Stage

KSPGMRESOrthog        13 1.0 4.7565e-04 4.4 1.44e+05 1.6 0.0e+00 0.0e+00 1.3e+01  0  5  0  0 10   0  5  0  0 10   496
KSPSetUp               2 1.0 3.6955e-05 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
KSPSolve               1 1.0 1.0792e-02 1.0 3.27e+06 2.4 2.8e+01 1.1e+03 3.7e+01  1100 41 13 28   1100 41 13 28   429
PCSetUp                2 1.0 9.3787e-03 3.0 7.17e+05 4.9 0.0e+00 0.0e+00 9.0e+00  1 19  0  0  7   1 19  0  0  7    92
PCSetUpOnBlocks        1 1.0 9.1140e-03 3.1 7.17e+05 4.9 0.0e+00 0.0e+00 7.0e+00  1 19  0  0  5   1 19  0  0  5    95
PCApply               15 1.0 6.8450e-04 1.9 1.74e+06 2.4 0.0e+00 0.0e+00 0.0e+00  0 53  0  0  0   0 53  0  0  0  3596
VecMDot               13 1.0 4.3201e-04 7.1 7.20e+04 1.6 0.0e+00 0.0e+00 1.3e+01  0  3  0  0 10   0  3  0  0 10   273
VecNorm               19 1.0 1.0467e-04 1.9 1.50e+04 1.6 0.0e+00 0.0e+00 1.9e+01  0  1  0  0 14   0  1  0  0 14   235
VecScale              14 1.0 1.6928e-05 1.1 5.54e+03 1.6 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0   536
VecCopy                5 1.0 6.1989e-06 1.4 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
VecSet                27 1.0 1.1921e-05 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
VecAXPY                3 1.0 1.0014e-05 1.2 2.38e+03 1.6 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0   388
VecMAXPY              14 1.0 2.3127e-05 1.4 8.24e+04 1.6 0.0e+00 0.0e+00 0.0e+00  0  3  0  0  0   0  3  0  0  0  5828
VecAssemblyBegin      18 1.0 2.9480e-0311.0 0.00e+00 0.0 4.0e+00 2.4e+03 4.5e+01  0  0  6  4 34   0  0  6  4 34     0
VecAssemblyEnd        18 1.0 2.7418e-05 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
VecScatterBegin       21 1.0 1.1706e-04 1.0 0.00e+00 0.0 4.2e+01 1.4e+03 0.0e+00  0  0 61 26  0   0  0 61 26  0     0
VecScatterEnd         21 1.0 6.2120e-03226.6 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
VecNormalize          14 1.0 8.5115e-05 1.2 1.66e+04 1.6 0.0e+00 0.0e+00 1.4e+01  0  1  0  0 11   0  1  0  0 11   320
MatMult               14 1.0 6.3713e-0323.5 6.38e+05 1.7 2.8e+01 1.1e+03 0.0e+00  0 22 41 13  0   0 22 41 13  0   161
MatSolve              15 1.0 5.5790e-04 2.4 1.74e+06 2.4 0.0e+00 0.0e+00 0.0e+00  0 53  0  0  0   0 53  0  0  0  4412
MatLUFactorNum         1 1.0 5.9414e-04 2.7 7.17e+05 4.9 0.0e+00 0.0e+00 0.0e+00  0 19  0  0  0   0 19  0  0  0  1455
MatILUFactorSym        1 1.0 8.3530e-03 3.3 0.00e+00 0.0 0.0e+00 0.0e+00 3.0e+00  1  0  0  0  2   1  0  0  0  2     0
MatAssemblyBegin       2 1.0 1.9193e-04 2.3 0.00e+00 0.0 3.0e+00 5.1e+04 4.0e+00  0  0  4 66  3   0  0  4 66  3     0
MatAssemblyEnd         2 1.0 8.3017e-04 1.2 0.00e+00 0.0 4.0e+00 2.8e+02 8.0e+00  0  0  6  0  6   0  0  6  0  6     0
MatGetRowIJ            1 1.0 2.8610e-06 1.5 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
MatGetOrdering         1 1.0 7.8917e-05 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 4.0e+00  0  0  0  0  3   0  0  0  0  3     0
MatZeroEntries         3 1.0 7.2002e-05 1.8 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
------------------------------------------------------------------------------------------------------------------------

Memory usage is given in bytes:

Object Type          Creations   Destructions     Memory  Descendants' Mem.
Reports information only for process 0.

--- Event Stage 0: Main Stage

       Krylov Solver     3              3        20592     0
      Preconditioner     3              3         2608     0
              Vector    38             38       169296     0
      Vector Scatter     5              5         5180     0
           Index Set    12             12        13760     0
   IS L to G Mapping     1              1          564     0
              Matrix     4              4      1034908     0
              Viewer     1              0            0     0
========================================================================================================================
Average time to get PetscTime(): 9.53674e-08
Average time for MPI_Barrier(): 2.19345e-06
Average time for zero size MPI_Send(): 1.34706e-05
#PETSc Option Table entries:
-ksp_right_pc
-log_summary
-pc_type bjacobi
-sub_pc_factor_levels 4
-sub_pc_factor_zeropivot 0
-sub_pc_type ilu
#End of PETSc Option Table entries
Compiled without FORTRAN kernels
Compiled with full precision matrices (default)
sizeof(short) 2 sizeof(int) 4 sizeof(long) 8 sizeof(void*) 8 sizeof(PetscScalar) 8 sizeof(PetscInt) 4
Configure run at: Thu Nov  8 11:21:02 2012
Configure options: --with-debugging=false --COPTFLAGS=-O3 --CXXOPTFLAGS=-O3 --FOPTFLAGS=-O3 --with-clanguage=C++ --with-shared-libraries=1 --with-mpi-dir=/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1 --with-mumps=true --download-mumps=1 --with-metis=true --download-metis=1 --with-parmetis=true --download-parmetis=1 --with-superlu=true --download-superlu=1 --with-superludir=true --download-superlu_dist=1 --with-blacs=true --download-blacs=1 --with-scalapack=true --download-scalapack=1 --with-hypre=true --download-hypre=1 --with-blas-lib="[/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_intel_lp64.so,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_sequential.so,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_core.so]" --with-lapack-lib="[/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_lapack95_lp64.a]"
-----------------------------------------
Libraries compiled on Thu Nov  8 11:21:02 2012 on daedalus.ices.utexas.edu 
Machine characteristics: Linux-2.6.32-279.1.1.el6.x86_64-x86_64-with-redhat-6.3-Carbon
Using PETSc directory: /opt/apps/ossw/libraries/petsc/petsc-3.3-p2
Using PETSc arch: intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt
-----------------------------------------

Using C compiler: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpicxx  -wd1572 -O3   -fPIC   ${COPTFLAGS} ${CFLAGS}
Using Fortran compiler: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpif90  -fPIC -O3   ${FOPTFLAGS} ${FFLAGS} 
-----------------------------------------

Using include paths: -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/include -I/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/include
-----------------------------------------

Using C linker: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpicxx
Using Fortran linker: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpif90
Using libraries: -Wl,-rpath,/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -L/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -lpetsc -lX11 -Wl,-rpath,/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -L/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -lcmumps -ldmumps -lsmumps -lzmumps -lmumps_common -lpord -lHYPRE -lpthread -lsuperlu_dist_3.0 -lparmetis -lmetis -lscalapack -lblacs -lsuperlu_4.3 -Wl,-rpath,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64 -L/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64 -lmkl_lapack95_lp64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -Wl,-rpath,/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/lib -L/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/lib -Wl,-rpath,/opt/apps/sysnet/intel/12.1/composer_xe_2011_sp1.7.256/compiler/lib/intel64 -L/opt/apps/sysnet/intel/12.1/composer_xe_2011_sp1.7.256/compiler/lib/intel64 -Wl,-rpath,/usr/lib/gcc/x86_64-redhat-linux/4.4.6 -L/usr/lib/gcc/x86_64-redhat-linux/4.4.6 -lmpichf90 -lifport -lifcore -lm -lm -lmpichcxx -ldl -lmpich -lopa -lmpl -lrt -lpthread -limf -lsvml -lipgo -ldecimal -lcilkrts -lstdc++ -lgcc_s -lirc -lirc_s -ldl 
-----------------------------------------


 ----------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                    |
| Num Processors: 2                                                                                                    |
| Time:           Tue Feb  5 13:54:13 2013                                                                             |
| OS:             Linux                                                                                                |
| HostName:       hbar.ices.utexas.edu                                                                                 |
| OS Release:     2.6.32-279.1.1.el6.x86_64                                                                            |
| OS Version:     #1 SMP Tue Jul 10 11:24:23 CDT 2012                                                                  |
| Machine:        x86_64                                                                                               |
| Username:       benkirk                                                                                              |
| Configuration:  ./configure  '--enable-everything'                                                                   |
|  '--prefix=/workspace/libmesh/install'                                                                               |
|  'CXX=mpicxx'                                                                                                        |
|  'CC=mpicc'                                                                                                          |
|  'F77=mpif77'                                                                                                        |
|  'FC=mpif90'                                                                                                         |
|  'PETSC_DIR=/opt/apps/ossw/libraries/petsc/petsc-3.3-p2'                                                             |
|  'PETSC_ARCH=intel-12.1-mkl-intel-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/intel-12.1/mpich2-1.4.1p1/mkl-intel-10.3.12.361'|
|  'VTK_DIR=/opt/apps/ossw/libraries/vtk/vtk-5.10.0/sl6/intel-12.1'                                                    |
 ----------------------------------------------------------------------------------------------------------------------
 ----------------------------------------------------------------------------------------------------------------
| libMesh Performance: Alive time=0.963037, Active time=0.887041                                                 |
 ----------------------------------------------------------------------------------------------------------------
| Event                              nCalls    Total Time  Avg Time    Total Time  Avg Time    % of Active Time  |
|                                              w/o Sub     w/o Sub     With Sub    With Sub    w/o S    With S   |
|----------------------------------------------------------------------------------------------------------------|
|                                                                                                                |
|                                                                                                                |
| DofMap                                                                                                         |
|   add_neighbors_to_send_list()     1         0.0060      0.006011    0.0158      0.015806    0.68     1.78     |
|   build_constraint_matrix()        124       0.0047      0.000038    0.0047      0.000038    0.53     0.53     |
|   build_sparsity()                 1         0.0159      0.015890    0.0286      0.028601    1.79     3.22     |
|   cnstrn_elem_mat_vec()            62        0.0083      0.000134    0.0083      0.000134    0.94     0.94     |
|   constrain_elem_vector()          62        0.0007      0.000012    0.0007      0.000012    0.08     0.08     |
|   create_dof_constraints()         1         0.0063      0.006289    0.0246      0.024643    0.71     2.78     |
|   distribute_dofs()                1         0.0023      0.002255    0.0059      0.005888    0.25     0.66     |
|   dof_indices()                    688       0.1020      0.000148    0.1020      0.000148    11.50    11.50    |
|   enforce_constraints_exactly()    3         0.0007      0.000244    0.0007      0.000244    0.08     0.08     |
|   prepare_send_list()              1         0.0002      0.000151    0.0002      0.000151    0.02     0.02     |
|   reinit()                         1         0.0034      0.003415    0.0034      0.003415    0.38     0.38     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          2         0.0009      0.000446    0.0196      0.009806    0.10     2.21     |
|                                                                                                                |
| ExodusII_IO                                                                                                    |
|   write_nodal_data()               1         0.0024      0.002415    0.0024      0.002415    0.27     0.27     |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        340       0.5813      0.001710    0.5813      0.001710    65.53    65.53    |
|   init_shape_functions()           157       0.0039      0.000025    0.0039      0.000025    0.44     0.44     |
|                                                                                                                |
| FEMSystem                                                                                                      |
|   assembly()                       1         0.0281      0.028111    0.0744      0.074414    3.17     8.39     |
|   assembly(get_residual)           1         0.0062      0.006215    0.0445      0.044471    0.70     5.01     |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             340       0.0058      0.000017    0.0058      0.000017    0.65     0.65     |
|   compute_face_map()               154       0.0019      0.000012    0.0019      0.000012    0.22     0.22     |
|   init_face_shape_functions()      2         0.0001      0.000034    0.0001      0.000034    0.01     0.01     |
|   init_reference_to_physical_map() 157       0.0045      0.000029    0.0045      0.000029    0.51     0.51     |
|                                                                                                                |
| Mesh                                                                                                           |
|   find_neighbors()                 1         0.0029      0.002855    0.0029      0.002941    0.32     0.33     |
|   renumber_nodes_and_elem()        2         0.0002      0.000085    0.0002      0.000085    0.02     0.02     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   compute_hilbert_indices()        2         0.0013      0.000633    0.0013      0.000633    0.14     0.14     |
|   find_global_indices()            2         0.0007      0.000330    0.0029      0.001448    0.07     0.33     |
|   parallel_sort()                  2         0.0006      0.000300    0.0007      0.000344    0.07     0.08     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         2         0.0014      0.000706    0.0236      0.011776    0.16     2.66     |
|                                                                                                                |
| MeshTools::Generation                                                                                          |
|   build_cube()                     1         0.0009      0.000908    0.0009      0.000908    0.10     0.10     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      1         0.0021      0.002103    0.0032      0.003228    0.24     0.36     |
|                                                                                                                |
| NewtonSolver                                                                                                   |
|   solve()                          1         0.0036      0.003586    0.1354      0.135409    0.40     15.27    |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      9         0.0002      0.000019    0.0002      0.000021    0.02     0.02     |
|   max(bool)                        1         0.0000      0.000003    0.0000      0.000003    0.00     0.00     |
|   max(scalar)                      150       0.0702      0.000468    0.0702      0.000468    7.91     7.91     |
|   max(vector)                      34        0.0002      0.000006    0.0004      0.000013    0.02     0.05     |
|   min(bool)                        175       0.0004      0.000002    0.0004      0.000002    0.04     0.04     |
|   min(scalar)                      143       0.0032      0.000022    0.0032      0.000022    0.36     0.36     |
|   min(vector)                      34        0.0003      0.000007    0.0005      0.000015    0.03     0.06     |
|   probe()                          12        0.0004      0.000032    0.0004      0.000032    0.04     0.04     |
|   receive()                        12        0.0001      0.000008    0.0005      0.000041    0.01     0.06     |
|   send()                           12        0.0000      0.000004    0.0000      0.000004    0.01     0.01     |
|   send_receive()                   16        0.0003      0.000016    0.0008      0.000052    0.03     0.09     |
|   sum()                            27        0.0002      0.000007    0.0003      0.000012    0.02     0.04     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           12        0.0000      0.000003    0.0000      0.000003    0.00     0.00     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         1         0.0004      0.000385    0.0005      0.000463    0.04     0.05     |
|   set_parent_processor_ids()       1         0.0001      0.000130    0.0001      0.000130    0.01     0.01     |
|                                                                                                                |
| PetscLinearSolver                                                                                              |
|   solve()                          1         0.0119      0.011924    0.0119      0.011924    1.34     1.34     |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            2754      0.8870                                          100.00            |
 ----------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example vector_fe_ex2:
*  mpirun -np 2 example-devel  -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************
Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:14 UTC

Hosted By:
SourceForge.net Logo