The source file exact_solution.C with comments:

        #include <math.h>
        
Mesh library includes
        #include "libmesh/libmesh_common.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
        
        
        
        
        /**
         *
         */
        Real exact_solution (const Real x,
        		     const Real y,
        		     const Real t)
        {
          const Real xo = 0.2;
          const Real yo = 0.2;
          const Real u  = 0.8;
          const Real v  = 0.8;
        
          const Real num =
            pow(x - u*t - xo, 2.) +
            pow(y - v*t - yo, 2.);
        
          const Real den =
            0.01*(4.*t + 1.);
        
          return exp(-num/den)/(4.*t + 1.);
        }



The source file transient_ex1.C with comments:

Transient Example 1 - Solving a Transient Linear System in Parallel



This example shows how a simple, linear transient system can be solved in parallel. The system is simple scalar convection-diffusion with a specified external velocity. The initial condition is given, and the solution is advanced in time with a standard Crank-Nicolson time-stepping strategy.

C++ include files that we need
        #include <iostream>
        #include <algorithm>
        #include <sstream>
        #include <math.h>
        
Basic include file needed for the mesh functionality.
        #include "libmesh/libmesh.h"
        #include "libmesh/mesh.h"
        #include "libmesh/mesh_refinement.h"
        #include "libmesh/gmv_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"
        
This example will solve a linear transient system, so we need to include the \p TransientLinearImplicitSystem definition.
        #include "libmesh/linear_implicit_system.h"
        #include "libmesh/transient_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 at each time step. Note that since the system is linear we technically do not need to assmeble the matrix at each time step, but we will anyway. In subsequent examples we will employ adaptive mesh refinement, and with a changing mesh it will be necessary to rebuild the system matrix.
        void assemble_cd (EquationSystems& es,
                          const std::string& system_name);
        
Function prototype. This function will initialize the system. Initialization functions are optional for systems. They allow you to specify the initial values of the solution. If an initialization function is not provided then the default (0) solution is provided.
        void init_cd (EquationSystems& es,
                      const std::string& system_name);
        
Exact solution function prototype. This gives the exact solution as a function of space and time. In this case the initial condition will be taken as the exact solution at time 0, as will the Dirichlet boundary conditions at time t.
        Real exact_solution (const Real x,
                             const Real y,
                             const Real t);
        
        Number exact_value (const Point& p,
                            const Parameters& parameters,
                            const std::string&,
                            const std::string&)
        {
          return exact_solution(p(0), p(1), parameters.get<Real> ("time"));
        }
        
        
        
We can now begin the main program. Note that this example will fail if you are using complex numbers since it was designed to be run only with real numbers.
        int main (int argc, char** argv)
        {
Initialize libMesh.
          LibMeshInit init (argc, argv);
        
This example requires Adaptive Mesh Refinement support - although it only refines uniformly, the refinement code used is the same underneath
        #ifndef LIBMESH_ENABLE_AMR
          libmesh_example_assert(false, "--enable-amr");
        #else
        
Skip this 2D example if libMesh was compiled as 1D-only.
          libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
        
Read the mesh from file. This is the coarse mesh that will be used in example 10 to demonstrate adaptive mesh refinement. Here we will simply read it in and uniformly refine it 5 times before we compute with it.

Create a mesh object, with dimension to be overridden later, distributed across the default MPI communicator.
          Mesh mesh(init.comm());
        
          mesh.read ("mesh.xda");
        
Create a MeshRefinement object to handle refinement of our mesh. This class handles all the details of mesh refinement and coarsening.
          MeshRefinement mesh_refinement (mesh);
        
Uniformly refine the mesh 5 times. This is the first time we use the mesh refinement capabilities of the library.
          mesh_refinement.uniformly_refine (5);
        
Print information about the mesh to the screen.
          mesh.print_info();
        
Create an equation systems object.
          EquationSystems equation_systems (mesh);
        
Add a transient system to the EquationSystems object named "Convection-Diffusion".
          TransientLinearImplicitSystem & system =
            equation_systems.add_system<TransientLinearImplicitSystem> ("Convection-Diffusion");
        
Adds the variable "u" to "Convection-Diffusion". "u" will be approximated using first-order approximation.
          system.add_variable ("u", FIRST);
        
Give the system a pointer to the matrix assembly and initialization functions.
          system.attach_assemble_function (assemble_cd);
          system.attach_init_function (init_cd);
        
Initialize the data structures for the equation system.
          equation_systems.init ();
        
Prints information about the system to the screen.
          equation_systems.print_info();
        
Write out the initial conditions.
          GMVIO(mesh).write_equation_systems ("out_000.gmv",
                                              equation_systems);
        
The Convection-Diffusion system requires that we specify the flow velocity. We will specify it as a RealVectorValue data type and then use the Parameters object to pass it to the assemble function.
          equation_systems.parameters.set<RealVectorValue>("velocity") =
            RealVectorValue (0.8, 0.8);
        
Solve the system "Convection-Diffusion". This will be done by looping over the specified time interval and calling the solve() member at each time step. This will assemble the system and call the linear solver.
          const Real dt = 0.025;
          system.time   = 0.;
        
          for (unsigned int t_step = 0; t_step < 50; t_step++)
            {
Incremenet the time counter, set the time and the time step size as parameters in the EquationSystem.
              system.time += dt;
        
              equation_systems.parameters.set<Real> ("time") = system.time;
              equation_systems.parameters.set<Real> ("dt")   = dt;
        
A pretty update message
              std::cout << " Solving time step ";
        
Since some compilers fail to offer full stream functionality, libMesh offers a string stream to work around this. Note that for other compilers, this is just a set of preprocessor macros and therefore should cost nothing (compared to a hand-coded string stream). We use additional curly braces here simply to enforce data locality.
              {
                std::ostringstream out;
        
                out << std::setw(2)
                    << std::right
                    << t_step
                    << ", time="
                    << std::fixed
                    << std::setw(6)
                    << std::setprecision(3)
                    << std::setfill('0')
                    << std::left
                    << system.time
                    <<  "...";
        
                std::cout << out.str() << std::endl;
              }
        
At this point we need to update the old solution vector. The old solution vector will be the current solution vector from the previous time step. We will do this by extracting the system from the \p EquationSystems object and using vector assignment. Since only \p TransientSystems (and systems derived from them) contain old solutions we need to specify the system type when we ask for it.
              *system.old_local_solution = *system.current_local_solution;
        
Assemble & solve the linear system
              equation_systems.get_system("Convection-Diffusion").solve();
        
Output evey 10 timesteps to file.
              if ( (t_step+1)%10 == 0)
                {
                  std::ostringstream file_name;
        
                  file_name << "out_"
                            << std::setw(3)
                            << std::setfill('0')
                            << std::right
                            << t_step+1
                            << ".gmv";
        
                  GMVIO(mesh).write_equation_systems (file_name.str(),
                                                      equation_systems);
                }
            }
        #endif // #ifdef LIBMESH_ENABLE_AMR
        
All done.
          return 0;
        }
        
We now define the function which provides the initialization routines for the "Convection-Diffusion" system. This handles things like setting initial conditions and boundary conditions.
        void init_cd (EquationSystems& es,
                      const std::string& system_name)
        {
It is a good idea to make sure we are initializing the proper system.
          libmesh_assert_equal_to (system_name, "Convection-Diffusion");
        
Get a reference to the Convection-Diffusion system object.
          TransientLinearImplicitSystem & system =
            es.get_system<TransientLinearImplicitSystem>("Convection-Diffusion");
        
Project initial conditions at time 0
          es.parameters.set<Real> ("time") = system.time = 0;
        
          system.project_solution(exact_value, NULL, es.parameters);
        }
        
        
        
Now we define the assemble function which will be used by the EquationSystems object at each timestep to assemble the linear system for solution.
        void assemble_cd (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, "Convection-Diffusion");
        
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.
          TransientLinearImplicitSystem & system =
            es.get_system<TransientLinearImplicitSystem> ("Convection-Diffusion");
        
Get a constant reference to 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<Point>& 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;
        
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;
        
Here we extract the velocity & parameters that we put in the EquationSystems object.
          const RealVectorValue velocity =
            es.parameters.get<RealVectorValue> ("velocity");
        
          const Real dt = es.parameters.get<Real>   ("dt");
        
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());
        
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++)
                {
Values to hold the old solution & its gradient.
                  Number   u_old = 0.;
                  Gradient grad_u_old;
        
Compute the old solution & its gradient.
                  for (unsigned int l=0; l<phi.size(); l++)
                    {
                      u_old      += phi[l][qp]*system.old_solution  (dof_indices[l]);
        
This will work, grad_u_old += dphi[l][qp]*system.old_solution (dof_indices[l]); but we can do it without creating a temporary like this:
                      grad_u_old.add_scaled (dphi[l][qp],system.old_solution (dof_indices[l]));
                    }
        
Now compute the element matrix and RHS contributions.
                  for (unsigned int i=0; i<phi.size(); i++)
                    {
The RHS contribution
                      Fe(i) += JxW[qp]*(
Mass matrix term
                                        u_old*phi[i][qp] +
                                        -.5*dt*(
Convection term (grad_u_old may be complex, so the order here is important!)
                                                (grad_u_old*velocity)*phi[i][qp] +
        
Diffusion term
                                                0.01*(grad_u_old*dphi[i][qp]))
                                        );
        
                      for (unsigned int j=0; j<phi.size(); j++)
                        {
The matrix contribution
                          Ke(i,j) += JxW[qp]*(
Mass-matrix
                                              phi[i][qp]*phi[j][qp] +
        
                                              .5*dt*(
Convection term
                                                     (velocity*dphi[j][qp])*phi[i][qp] +
        
Diffusion term
                                                     0.01*(dphi[i][qp]*dphi[j][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++)
                        {
                          const Number value = exact_solution (qface_points[qp](0),
                                                               qface_points[qp](1),
                                                               system.time);
        
RHS contribution
                          for (unsigned int i=0; i<psi.size(); i++)
                            Fe(i) += penalty*JxW_face[qp]*value*psi[i][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];
                        }
                    }
              }
        
If this assembly program were to be used on an adaptive mesh, we would have to apply any hanging node constraint equations
              dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
        
The element matrix and right-hand-side are now built for this element. Add them to the global matrix and right-hand-side vector. The \p SparseMatrix::add_matrix() and \p NumericVector::add_vector() members do this for us.
              system.matrix->add_matrix (Ke, dof_indices);
              system.rhs->add_vector    (Fe, dof_indices);
            }
        
That concludes the system matrix assembly routine.
        #endif // #ifdef LIBMESH_ENABLE_AMR
        }



The source file exact_solution.C without comments:

 
  #include <math.h>
  
  #include "libmesh/libmesh_common.h"
  
  using namespace libMesh;
  
  
  
  
  
  /**
   *
   */
  Real exact_solution (const Real x,
  		     const Real y,
  		     const Real t)
  {
    const Real xo = 0.2;
    const Real yo = 0.2;
    const Real u  = 0.8;
    const Real v  = 0.8;
  
    const Real num =
      pow(x - u*t - xo, 2.) +
      pow(y - v*t - yo, 2.);
  
    const Real den =
      0.01*(4.*t + 1.);
  
    return exp(-num/den)/(4.*t + 1.);
  }



The source file transient_ex1.C without comments:

 
  
  #include <iostream>
  #include <algorithm>
  #include <sstream>
  #include <math.h>
  
  #include "libmesh/libmesh.h"
  #include "libmesh/mesh.h"
  #include "libmesh/mesh_refinement.h"
  #include "libmesh/gmv_io.h"
  #include "libmesh/equation_systems.h"
  #include "libmesh/fe.h"
  #include "libmesh/quadrature_gauss.h"
  #include "libmesh/dof_map.h"
  #include "libmesh/sparse_matrix.h"
  #include "libmesh/numeric_vector.h"
  #include "libmesh/dense_matrix.h"
  #include "libmesh/dense_vector.h"
  
  #include "libmesh/linear_implicit_system.h"
  #include "libmesh/transient_system.h"
  #include "libmesh/vector_value.h"
  
  #include "libmesh/elem.h"
  
  using namespace libMesh;
  
  void assemble_cd (EquationSystems& es,
                    const std::string& system_name);
  
  void init_cd (EquationSystems& es,
                const std::string& system_name);
  
  Real exact_solution (const Real x,
                       const Real y,
                       const Real t);
  
  Number exact_value (const Point& p,
                      const Parameters& parameters,
                      const std::string&,
                      const std::string&)
  {
    return exact_solution(p(0), p(1), parameters.get<Real> ("time"));
  }
  
  
  
  int main (int argc, char** argv)
  {
    LibMeshInit init (argc, argv);
  
  #ifndef LIBMESH_ENABLE_AMR
    libmesh_example_assert(false, "--enable-amr");
  #else
  
    libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
  
    Mesh mesh(init.comm());
  
    mesh.read ("mesh.xda");
  
    MeshRefinement mesh_refinement (mesh);
  
    mesh_refinement.uniformly_refine (5);
  
    mesh.print_info();
  
    EquationSystems equation_systems (mesh);
  
    TransientLinearImplicitSystem & system =
      equation_systems.add_system<TransientLinearImplicitSystem> ("Convection-Diffusion");
  
    system.add_variable ("u", FIRST);
  
    system.attach_assemble_function (assemble_cd);
    system.attach_init_function (init_cd);
  
    equation_systems.init ();
  
    equation_systems.print_info();
  
    GMVIO(mesh).write_equation_systems ("out_000.gmv",
                                        equation_systems);
  
    equation_systems.parameters.set<RealVectorValue>("velocity") =
      RealVectorValue (0.8, 0.8);
  
    const Real dt = 0.025;
    system.time   = 0.;
  
    for (unsigned int t_step = 0; t_step < 50; t_step++)
      {
        system.time += dt;
  
        equation_systems.parameters.set<Real> ("time") = system.time;
        equation_systems.parameters.set<Real> ("dt")   = dt;
  
        std::cout << " Solving time step ";
  
        {
          std::ostringstream out;
  
          out << std::setw(2)
              << std::right
              << t_step
              << ", time="
              << std::fixed
              << std::setw(6)
              << std::setprecision(3)
              << std::setfill('0')
              << std::left
              << system.time
              <<  "...";
  
          std::cout << out.str() << std::endl;
        }
  
        *system.old_local_solution = *system.current_local_solution;
  
        equation_systems.get_system("Convection-Diffusion").solve();
  
        if ( (t_step+1)%10 == 0)
          {
            std::ostringstream file_name;
  
            file_name << "out_"
                      << std::setw(3)
                      << std::setfill('0')
                      << std::right
                      << t_step+1
                      << ".gmv";
  
            GMVIO(mesh).write_equation_systems (file_name.str(),
                                                equation_systems);
          }
      }
  #endif // #ifdef LIBMESH_ENABLE_AMR
  
    return 0;
  }
  
  void init_cd (EquationSystems& es,
                const std::string& system_name)
  {
    libmesh_assert_equal_to (system_name, "Convection-Diffusion");
  
    TransientLinearImplicitSystem & system =
      es.get_system<TransientLinearImplicitSystem>("Convection-Diffusion");
  
    es.parameters.set<Real> ("time") = system.time = 0;
  
    system.project_solution(exact_value, NULL, es.parameters);
  }
  
  
  
  void assemble_cd (EquationSystems& es,
                    const std::string& system_name)
  {
  #ifdef LIBMESH_ENABLE_AMR
    libmesh_assert_equal_to (system_name, "Convection-Diffusion");
  
    const MeshBase& mesh = es.get_mesh();
  
    const unsigned int dim = mesh.mesh_dimension();
  
    TransientLinearImplicitSystem & system =
      es.get_system<TransientLinearImplicitSystem> ("Convection-Diffusion");
  
    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 std::vector<Point>& qface_points = fe_face->get_xyz();
  
    const DofMap& dof_map = system.get_dof_map();
  
    DenseMatrix<Number> Ke;
    DenseVector<Number> Fe;
  
    std::vector<dof_id_type> dof_indices;
  
    const RealVectorValue velocity =
      es.parameters.get<RealVectorValue> ("velocity");
  
    const Real dt = es.parameters.get<Real>   ("dt");
  
    MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
    const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
  
    for ( ; el != end_el; ++el)
      {
        const Elem* elem = *el;
  
        dof_map.dof_indices (elem, dof_indices);
  
        fe->reinit (elem);
  
        Ke.resize (dof_indices.size(),
                   dof_indices.size());
  
        Fe.resize (dof_indices.size());
  
        for (unsigned int qp=0; qp<qrule.n_points(); qp++)
          {
            Number   u_old = 0.;
            Gradient grad_u_old;
  
            for (unsigned int l=0; l<phi.size(); l++)
              {
                u_old      += phi[l][qp]*system.old_solution  (dof_indices[l]);
  
                grad_u_old.add_scaled (dphi[l][qp],system.old_solution (dof_indices[l]));
              }
  
            for (unsigned int i=0; i<phi.size(); i++)
              {
                Fe(i) += JxW[qp]*(
                                  u_old*phi[i][qp] +
                                  -.5*dt*(
                                          (grad_u_old*velocity)*phi[i][qp] +
  
                                          0.01*(grad_u_old*dphi[i][qp]))
                                  );
  
                for (unsigned int j=0; j<phi.size(); j++)
                  {
                    Ke(i,j) += JxW[qp]*(
                                        phi[i][qp]*phi[j][qp] +
  
                                        .5*dt*(
                                               (velocity*dphi[j][qp])*phi[i][qp] +
  
                                               0.01*(dphi[i][qp]*dphi[j][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++)
                  {
                    const Number value = exact_solution (qface_points[qp](0),
                                                         qface_points[qp](1),
                                                         system.time);
  
                    for (unsigned int i=0; i<psi.size(); i++)
                      Fe(i) += penalty*JxW_face[qp]*value*psi[i][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];
                  }
              }
        }
  
        dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
  
        system.matrix->add_matrix (Ke, dof_indices);
        system.rhs->add_vector    (Fe, dof_indices);
      }
  
  #endif // #ifdef LIBMESH_ENABLE_AMR
  }



The console output of the program:

make[4]: Entering directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/transient/transient_ex1'
***************************************************************
* Running Example transient_ex1:
*  mpirun -np 4 example-devel  -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
 
 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=6273
    n_local_nodes()=1789
  n_elem()=13650
    n_local_elem()=3277
    n_active_elem()=10240
  n_subdomains()=1
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System #0, "Convection-Diffusion"
    Type "TransientLinearImplicit"
    Variables="u" 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="FIRST", "THIRD" 
    n_dofs()=6273
    n_local_dofs()=1789
    n_constrained_dofs()=0
    n_local_constrained_dofs()=0
    n_vectors()=3
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 7.53164
      Average Off-Processor Bandwidth <= 0.127212
      Maximum  On-Processor Bandwidth <= 11
      Maximum Off-Processor Bandwidth <= 5
    DofMap Constraints
      Number of DoF Constraints = 0
      Number of Node Constraints = 0

 Solving time step  0, time=0.0250...
 Solving time step  1, time=0.0500...
 Solving time step  2, time=0.0750...
 Solving time step  3, time=0.1000...
 Solving time step  4, time=0.1250...
 Solving time step  5, time=0.1500...
 Solving time step  6, time=0.1750...
 Solving time step  7, time=0.2000...
 Solving time step  8, time=0.2250...
 Solving time step  9, time=0.2500...
 Solving time step 10, time=0.2750...
 Solving time step 11, time=0.3000...
 Solving time step 12, time=0.3250...
 Solving time step 13, time=0.3500...
 Solving time step 14, time=0.3750...
 Solving time step 15, time=0.4000...
 Solving time step 16, time=0.4250...
 Solving time step 17, time=0.4500...
 Solving time step 18, time=0.4750...
 Solving time step 19, time=0.5000...
 Solving time step 20, time=0.5250...
 Solving time step 21, time=0.5500...
 Solving time step 22, time=0.5750...
 Solving time step 23, time=0.6000...
 Solving time step 24, time=0.6250...
 Solving time step 25, time=0.6500...
 Solving time step 26, time=0.6750...
 Solving time step 27, time=0.7000...
 Solving time step 28, time=0.7250...
 Solving time step 29, time=0.7500...
 Solving time step 30, time=0.7750...
 Solving time step 31, time=0.8000...
 Solving time step 32, time=0.8250...
 Solving time step 33, time=0.8500...
 Solving time step 34, time=0.8750...
 Solving time step 35, time=0.9000...
 Solving time step 36, time=0.9250...
 Solving time step 37, time=0.9500...
 Solving time step 38, time=0.9750...
 Solving time step 39, time=1.0000...
 Solving time step 40, time=1.0250...
 Solving time step 41, time=1.0500...
 Solving time step 42, time=1.0750...
 Solving time step 43, time=1.1000...
 Solving time step 44, time=1.1250...
 Solving time step 45, time=1.1500...
 Solving time step 46, time=1.1750...
 Solving time step 47, time=1.2000...
 Solving time step 48, time=1.2250...
 Solving time step 49, time=1.2500...

 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:57:44 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=3.94001, Active time=3.80371                                                   |
 ----------------------------------------------------------------------------------------------------------------
| 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.0118      0.011765    0.0132      0.013152    0.31     0.35     |
|   build_sparsity()                 1         0.0080      0.007973    0.0200      0.019989    0.21     0.53     |
|   create_dof_constraints()         1         0.0134      0.013400    0.0134      0.013400    0.35     0.35     |
|   distribute_dofs()                1         0.0241      0.024060    0.0638      0.063845    0.63     1.68     |
|   dof_indices()                    139514    0.5224      0.000004    0.5224      0.000004    13.73    13.73    |
|   prepare_send_list()              1         0.0000      0.000021    0.0000      0.000021    0.00     0.00     |
|   reinit()                         1         0.0392      0.039237    0.0392      0.039237    1.03     1.03     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          6         0.0354      0.005902    0.1033      0.017222    0.93     2.72     |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        123800    0.2580      0.000002    0.2580      0.000002    6.78     6.78     |
|   init_shape_functions()           3900      0.0053      0.000001    0.0053      0.000001    0.14     0.14     |
|   inverse_map()                    7600      0.0217      0.000003    0.0217      0.000003    0.57     0.57     |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             123800    0.2322      0.000002    0.2322      0.000002    6.11     6.11     |
|   compute_face_map()               3800      0.0179      0.000005    0.0403      0.000011    0.47     1.06     |
|   init_face_shape_functions()      100       0.0003      0.000003    0.0003      0.000003    0.01     0.01     |
|   init_reference_to_physical_map() 3900      0.0126      0.000003    0.0126      0.000003    0.33     0.33     |
|                                                                                                                |
| GMVIO                                                                                                          |
|   write_nodal_data()               6         0.2407      0.040114    0.2411      0.040185    6.33     6.34     |
|                                                                                                                |
| LocationMap                                                                                                    |
|   find()                           32736     0.0391      0.000001    0.0391      0.000001    1.03     1.03     |
|   init()                           5         0.0009      0.000171    0.0009      0.000171    0.02     0.02     |
|                                                                                                                |
| Mesh                                                                                                           |
|   find_neighbors()                 2         0.0484      0.024205    0.0487      0.024327    1.27     1.28     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   broadcast()                      1         0.0003      0.000256    0.0005      0.000544    0.01     0.01     |
|   compute_hilbert_indices()        3         0.0424      0.014144    0.0424      0.014144    1.12     1.12     |
|   find_global_indices()            3         0.0065      0.002163    0.0506      0.016862    0.17     1.33     |
|   parallel_sort()                  3         0.0009      0.000296    0.0010      0.000342    0.02     0.03     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         6         0.0002      0.000036    0.3450      0.057504    0.01     9.07     |
|                                                                                                                |
| MeshRefinement                                                                                                 |
|   _refine_elements()               5         0.0576      0.011520    0.1587      0.031734    1.51     4.17     |
|   add_point()                      32736     0.0548      0.000002    0.0975      0.000003    1.44     2.56     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      2         0.1815      0.090736    0.2316      0.115795    4.77     6.09     |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      11        0.0002      0.000015    0.0002      0.000019    0.00     0.01     |
|   broadcast()                      27        0.0001      0.000004    0.0001      0.000003    0.00     0.00     |
|   max(bool)                        6         0.0001      0.000012    0.0001      0.000012    0.00     0.00     |
|   max(scalar)                      272       0.0012      0.000004    0.0012      0.000004    0.03     0.03     |
|   max(vector)                      61        0.0004      0.000006    0.0011      0.000017    0.01     0.03     |
|   min(bool)                        322       0.0011      0.000003    0.0011      0.000003    0.03     0.03     |
|   min(scalar)                      264       0.0149      0.000056    0.0149      0.000056    0.39     0.39     |
|   min(vector)                      61        0.0005      0.000008    0.0012      0.000019    0.01     0.03     |
|   probe()                          48        0.0003      0.000007    0.0003      0.000007    0.01     0.01     |
|   receive()                        48        0.0003      0.000006    0.0006      0.000013    0.01     0.02     |
|   send()                           48        0.0002      0.000005    0.0002      0.000005    0.01     0.01     |
|   send_receive()                   54        0.0003      0.000006    0.0013      0.000024    0.01     0.03     |
|   sum()                            39        0.0018      0.000045    0.0150      0.000385    0.05     0.40     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           48        0.0001      0.000001    0.0001      0.000001    0.00     0.00     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         2         0.0099      0.004946    0.0102      0.005124    0.26     0.27     |
|   set_parent_processor_ids()       2         0.0067      0.003362    0.0067      0.003362    0.18     0.18     |
|                                                                                                                |
| PetscLinearSolver                                                                                              |
|   solve()                          50        0.7957      0.015914    0.7957      0.015914    20.92    20.92    |
|                                                                                                                |
| System                                                                                                         |
|   assemble()                       50        1.0811      0.021622    2.1137      0.042274    28.42    55.57    |
|   project_vector()                 1         0.0134      0.013423    0.0250      0.024956    0.35     0.66     |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            473348    3.8037                                          100.00            |
 ----------------------------------------------------------------------------------------------------------------

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

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

Hosted By:
SourceForge.net Logo