The source file transient_ex2.C with comments:

Transient Example 2 - The Newmark System and the Wave Equation



This is the eighth example program. It builds on the previous example programs. It introduces the NewmarkSystem class. In this example the wave equation is solved using the time integration scheme provided by the NewmarkSystem class.

This example comes with a cylindrical mesh given in the universal file pipe-mesh.unv. The mesh contains HEX8 and PRISM6 elements.

C++ include files that we need
        #include <iostream>
        #include <fstream>
        #include <algorithm>
        #include <stdio.h>
        #include <math.h>
        
Basic include file needed for the mesh functionality.
        #include "libmesh/libmesh.h"
        #include "libmesh/serial_mesh.h"
        #include "libmesh/gmv_io.h"
        #include "libmesh/vtk_io.h"
        #include "libmesh/newmark_system.h"
        #include "libmesh/equation_systems.h"
        
Define the Finite Element object.
        #include "libmesh/fe.h"
        
Define Gauss quadrature rules.
        #include "libmesh/quadrature_gauss.h"
        
Define useful datatypes for finite element
        #include "libmesh/dense_matrix.h"
        #include "libmesh/dense_vector.h"
        
Define matrix and vector data types for the global equation system. These are base classes, from which specific implementations, like the PETSc or LASPACK implementations, are derived.
        #include "libmesh/sparse_matrix.h"
        #include "libmesh/numeric_vector.h"
        
Define the DofMap, which handles degree of freedom indexing.
        #include "libmesh/dof_map.h"
        
The definition of a vertex associated with a Mesh.
        #include "libmesh/node.h"
        
The definition of a geometric element
        #include "libmesh/elem.h"
        
Defines the MeshData class, which allows you to store data about the mesh when reading in files, etc.
        #include "libmesh/mesh_data.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
Function prototype. This is the function that will assemble the linear system for our problem, governed by the linear wave equation.
        void assemble_wave(EquationSystems& es,
                           const std::string& system_name);
        
        
Function Prototype. This function will be used to apply the initial conditions.
        void apply_initial(EquationSystems& es,
                           const std::string& system_name);
        
Function Prototype. This function imposes Dirichlet Boundary conditions via the penalty method after the system is assembled.
        void fill_dirichlet_bc(EquationSystems& es,
                               const std::string& system_name);
        
The main program
        int main (int argc, char** argv)
        {
Initialize libraries, like in example 2.
          LibMeshInit init (argc, argv);
        
Check for proper usage.
          if (argc < 2)
            {
              if (libMesh::processor_id() == 0)
                std::cerr << "Usage: " << argv[0] << " [meshfile]"
                          << std::endl;
        
              libmesh_error();
            }
        
Tell the user what we are doing.
          else
            {
              std::cout << "Running " << argv[0];
        
              for (int i=1; i<argc; i++)
                std::cout << " " << argv[i];
        
              std::cout << std::endl << std::endl;
        
            }
        
LasPack solvers don't work so well for this example (not sure why), and Trilinos matrices don't work at all. Print a warning to the user if PETSc is not in use.
          if (libMesh::default_solver_package() == LASPACK_SOLVERS)
            {
              std::cout << "WARNING! It appears you are using the\n"
                        << "LasPack solvers.  This example may not converge\n"
                        << "using LasPack, but should work OK with PETSc.\n"
                        << "http://www.mcs.anl.gov/petsc/\n"
                        << std::endl;
            }
          else if (libMesh::default_solver_package() == TRILINOS_SOLVERS)
            {
              std::cout << "WARNING! It appears you are using the\n"
                        << "Trilinos solvers.  The current libMesh Epetra\n"
                        << "interface does not allow sparse matrix addition,\n"
                        << "as is needed in this problem.  We recommend\n"
                        << "using PETSc: http://www.mcs.anl.gov/petsc/\n"
                        << std::endl;
              return 0;
            }
        
Get the name of the mesh file from the command line.
          std::string mesh_file = argv[1];
          std::cout << "Mesh file is: " << mesh_file << std::endl;
        
Skip this 3D example if libMesh was compiled as 1D or 2D-only.
          libmesh_example_assert(3 <= LIBMESH_DIM, "3D support");
        
Create a mesh. This example directly references all mesh nodes and is incompatible with ParallelMesh use.

Create a SerialMesh object, with dimension to be overridden later, distributed across the default MPI communicator.
          SerialMesh mesh(init.comm());
          MeshData mesh_data(mesh);
        
Read the meshfile specified in the command line or use the internal mesh generator to create a uniform grid on an elongated cube.
          mesh.read(mesh_file, &mesh_data);
        
mesh.build_cube (10, 10, 40, -1., 1., -1., 1., 0., 4., HEX8);

Print information about the mesh to the screen.
          mesh.print_info();
        
The node that should be monitored.
          const unsigned int result_node = 274;
        
        
Time stepping issues

Note that the total current time is stored as a parameter in the \pEquationSystems object.

the time step size
          const Real delta_t = .0000625;
        
The number of time steps.
          unsigned int n_time_steps = 300;
        
Create an equation systems object.
          EquationSystems equation_systems (mesh);
        
Declare the system and its variables. Create a NewmarkSystem named "Wave"
          equation_systems.add_system<NewmarkSystem> ("Wave");
        
Use a handy reference to this system
          NewmarkSystem & t_system = equation_systems.get_system<NewmarkSystem> ("Wave");
        
Add the variable "p" to "Wave". "p" will be approximated using first-order approximation.
          t_system.add_variable("p", FIRST);
        
Give the system a pointer to the matrix assembly function and the initial condition function defined below.
          t_system.attach_assemble_function  (assemble_wave);
          t_system.attach_init_function      (apply_initial);
        
Set the time step size, and optionally the Newmark parameters, so that \p NewmarkSystem can compute integration constants. Here we simply use pass only the time step and use default values for \p alpha=.25 and \p delta=.5.
          t_system.set_newmark_parameters(delta_t);
        
Set the speed of sound and fluid density as \p EquationSystems parameter, so that \p assemble_wave() can access it.
          equation_systems.parameters.set<Real>("speed")          = 1000.;
          equation_systems.parameters.set<Real>("fluid density")  = 1000.;
        
Start time integration from t=0
          t_system.time = 0.;
        
Initialize the data structures for the equation system.
          equation_systems.init();
        
Prints information about the system to the screen.
          equation_systems.print_info();
        
A file to store the results at certain nodes.
          std::ofstream res_out("pressure_node.res");
        
get the dof_numbers for the nodes that should be monitored.
          const unsigned int res_node_no = result_node;
          const Node& res_node = mesh.node(res_node_no-1);
          unsigned int dof_no = res_node.dof_number(0,0,0);
        
Assemble the time independent system matrices and rhs. This function will also compute the effective system matrix K~=K+a_0*M+a_1*C and apply user specified initial conditions.
          t_system.assemble();
        
Now solve for each time step. For convenience, use a local buffer of the current time. But once this time is updated, also update the \p EquationSystems parameter Start with t_time = 0 and write a short header to the nodal result file
          res_out << "# pressure at node " << res_node_no << "\n"
                  << "# time\tpressure\n"
                  << t_system.time << "\t" << 0 << std::endl;
        
        
          for (unsigned int time_step=0; time_step<n_time_steps; time_step++)
            {
Update the time. Both here and in the \p EquationSystems object
              t_system.time += delta_t;
        
Update the rhs.
              t_system.update_rhs();
        
Impose essential boundary conditions. Not that since the matrix is only assembled once, the penalty parameter should be added to the matrix only in the first time step. The applied boundary conditions may be time-dependent and hence the rhs vector is considered in each time step.
              if (time_step == 0)
                {
The local function \p fill_dirichlet_bc() may also set Dirichlet boundary conditions for the matrix. When you set the flag as shown below, the flag will return true. If you want it to return false, simply do not set it.
                  equation_systems.parameters.set<bool>("Newmark set BC for Matrix") = true;
        
                  fill_dirichlet_bc(equation_systems, "Wave");
        
unset the flag, so that it returns false
                  equation_systems.parameters.set<bool>("Newmark set BC for Matrix") = false;
                }
              else
                fill_dirichlet_bc(equation_systems, "Wave");
        
Solve the system "Wave".
              t_system.solve();
        
After solving the system, write the solution to a GMV-formatted plot file. Do only for a few time steps.
              if (time_step == 30 || time_step == 60 ||
                  time_step == 90 || time_step == 120 )
                {
                  char buf[14];
        
        		  if (!libMesh::on_command_line("--vtk")){
        			  sprintf (buf, "out.%03d.gmv", time_step);
        
        	          GMVIO(mesh).write_equation_systems (buf,equation_systems);
        
        		  }else{
        #ifdef LIBMESH_HAVE_VTK
VTK viewers are generally not happy with two dots in a filename
                                  sprintf (buf, "out_%03d.exd", time_step);
        
        	          VTKIO(mesh).write_equation_systems (buf,equation_systems);
        #endif // #ifdef LIBMESH_HAVE_VTK
        		  }
                }
        
Update the p, v and a.
              t_system.update_u_v_a();
        
dof_no may not be local in parallel runs, so we may need a global displacement vector
              NumericVector<Number> &displacement
                = t_system.get_vector("displacement");
              std::vector<Number> global_displacement(displacement.size());
              displacement.localize(global_displacement);
        
Write nodal results to file. The results can then be viewed with e.g. gnuplot (run gnuplot and type 'plot "pressure_node.res" with lines' in the command line)
              res_out << t_system.time << "\t"
                      << global_displacement[dof_no]
                      << std::endl;
            }
        
All done.
          return 0;
        }
        
This function assembles the system matrix and right-hand-side for our wave equation.
        void assemble_wave(EquationSystems& es,
                           const std::string& system_name)
        {
It is a good idea to make sure we are assembling the proper system.
          libmesh_assert_equal_to (system_name, "Wave");
        
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();
        
Copy the speed of sound and fluid density to a local variable.
          const Real speed = es.parameters.get<Real>("speed");
          const Real rho   = es.parameters.get<Real>("fluid density");
        
Get a reference to our system, as before.
          NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);
        
Get a constant reference to the Finite Element type for the first (and only) variable in the system.
          FEType fe_type = t_system.get_dof_map().variable_type(0);
        
In here, we will add the element matrices to the @e additional matrices "stiffness_mass" and "damping" and the additional vector "force", not to the members "matrix" and "rhs". Therefore, get writable references to them.
          SparseMatrix<Number>&   stiffness = t_system.get_matrix("stiffness");
          SparseMatrix<Number>&   damping   = t_system.get_matrix("damping");
          SparseMatrix<Number>&   mass      = t_system.get_matrix("mass");
        
          NumericVector<Number>&  force     = t_system.get_vector("force");
        
Some solver packages (PETSc) are especially picky about allocating sparsity structure and truly assigning values to this structure. Namely, matrix additions, as performed later, exhibit acceptable performance only for identical sparsity structures. Therefore, explicitly zero the values in the collective matrix, so that matrix additions encounter identical sparsity structures.
          SparseMatrix<Number>&  matrix     = *t_system.matrix;
          DenseMatrix<Number>    zero_matrix;
        
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));
        
A 2nd order Gauss quadrature rule for numerical integration.
          QGauss qrule (dim, SECOND);
        
Tell the finite element object to use our quadrature rule.
          fe->attach_quadrature_rule (&qrule);
        
The element Jacobian * quadrature weight at each integration point.
          const std::vector<Real>& JxW = fe->get_JxW();
        
The element shape functions evaluated at the quadrature points.
          const std::vector<std::vector<Real> >& phi = fe->get_phi();
        
The element shape function gradients evaluated at the quadrature points.
          const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
        
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.
          const DofMap& dof_map = t_system.get_dof_map();
        
The element mass, damping and stiffness matrices and the element contribution to the rhs.
          DenseMatrix<Number>   Ke, Ce, Me;
          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;
        
Now we will loop over all the elements in the mesh. We will compute the element matrix and right-hand-side contribution.
          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 matrices and rhs 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 HEX8 and now have a PRISM6).
              {
                const unsigned int n_dof_indices = dof_indices.size();
        
                Ke.resize          (n_dof_indices, n_dof_indices);
                Ce.resize          (n_dof_indices, n_dof_indices);
                Me.resize          (n_dof_indices, n_dof_indices);
                zero_matrix.resize (n_dof_indices, n_dof_indices);
                Fe.resize          (n_dof_indices);
              }
        
Now loop over the quadrature points. This handles the numeric integration.
              for (unsigned int qp=0; qp<qrule.n_points(); qp++)
                {
Now we will build the element matrix. This involves a double loop to integrate the test funcions (i) against the trial functions (j).
                  for (unsigned int i=0; i<phi.size(); i++)
                    for (unsigned int j=0; j<phi.size(); j++)
                      {
                        Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
                        Me(i,j) += JxW[qp]*phi[i][qp]*phi[j][qp]
                                   *1./(speed*speed);
                      } // end of the matrix summation loop
                } // end of quadrature point loop
        
Now compute the contribution to the element matrix and the right-hand-side vector if the current element lies on the boundary.
              {
In this example no natural boundary conditions will be considered. The code is left here so it can easily be extended.

don't do this for any side
                for (unsigned int side=0; side<elem->n_sides(); side++)
                  if (!true)
if (elem->neighbor(side) == NULL)
                    {
Declare a special finite element object for boundary integration.
                      AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
        
Boundary integration requires one quadraure rule, with dimensionality one less than the dimensionality of the element.
                      QGauss qface(dim-1, SECOND);
        
Tell the finte element object to use our quadrature rule.
                      fe_face->attach_quadrature_rule (&qface);
        
The value of the shape functions at the quadrature points.
                      const std::vector<std::vector<Real> >&  phi_face = fe_face->get_phi();
        
The Jacobian * Quadrature Weight at the quadrature points on the face.
                      const std::vector<Real>& JxW_face = fe_face->get_JxW();
        
Compute the shape function values on the element face.
                      fe_face->reinit(elem, side);
        
Here we consider a normal acceleration acc_n=1 applied to the whole boundary of our mesh.
                      const Real acc_n_value = 1.0;
        
Loop over the face quadrature points for integration.
                      for (unsigned int qp=0; qp<qface.n_points(); qp++)
                        {
Right-hand-side contribution due to prescribed normal acceleration.
                          for (unsigned int i=0; i<phi_face.size(); i++)
                            {
                              Fe(i) += acc_n_value*rho
                                *phi_face[i][qp]*JxW_face[qp];
                            }
                        } // end face quadrature point loop
                    } // end if (elem->neighbor(side) == NULL)
        
In this example the Dirichlet boundary conditions will be imposed via panalty method after the system is assembled.

              } // end boundary condition section
        
If this assembly program were to be used on an adaptive mesh, we would have to apply any hanging node constraint equations by uncommenting the following lines: std::vector dof_indicesC = dof_indices; std::vector dof_indicesM = dof_indices; dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices); dof_map.constrain_element_matrix (Ce, dof_indicesC); dof_map.constrain_element_matrix (Me, dof_indicesM);

Finally, simply add the contributions to the additional matrices and vector.
              stiffness.add_matrix (Ke, dof_indices);
              damping.add_matrix   (Ce, dof_indices);
              mass.add_matrix      (Me, dof_indices);
        
              force.add_vector     (Fe, dof_indices);
        
For the overall matrix, explicitly zero the entries where we added values in the other ones, so that we have identical sparsity footprints.
              matrix.add_matrix(zero_matrix, dof_indices);
        
            } // end of element loop
        
All done!
          return;
        }
        
This function applies the initial conditions
        void apply_initial(EquationSystems& es,
                           const std::string& system_name)
        {
Get a reference to our system, as before
          NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);
        
Numeric vectors for the pressure, velocity and acceleration values.
          NumericVector<Number>&  pres_vec       = t_system.get_vector("displacement");
          NumericVector<Number>&  vel_vec        = t_system.get_vector("velocity");
          NumericVector<Number>&  acc_vec        = t_system.get_vector("acceleration");
        
Assume our fluid to be at rest, which would also be the default conditions in class NewmarkSystem, but let us do it explicetly here.
          pres_vec.zero();
          vel_vec.zero();
          acc_vec.zero();
        }
        
This function applies the Dirichlet boundary conditions
        void fill_dirichlet_bc(EquationSystems& es,
                               const std::string& system_name)
        {
It is a good idea to make sure we are assembling the proper system.
          libmesh_assert_equal_to (system_name, "Wave");
        
Get a reference to our system, as before.
          NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);
        
Get writable references to the overall matrix and vector.
          SparseMatrix<Number>&  matrix = *t_system.matrix;
          NumericVector<Number>& rhs    = *t_system.rhs;
        
Get a constant reference to the mesh object.
          const MeshBase& mesh = es.get_mesh();
        
Get \p libMesh's \f$ \pi \f$
          const Real pi = libMesh::pi;
        
Ask the \p EquationSystems flag whether we should do this also for the matrix
          const bool do_for_matrix =
            es.parameters.get<bool>("Newmark set BC for Matrix");
        
Number of nodes in the mesh.
          unsigned int n_nodes = mesh.n_nodes();
        
          for (unsigned int n_cnt=0; n_cnt<n_nodes; n_cnt++)
            {
Get a reference to the current node.
              const Node& curr_node = mesh.node(n_cnt);
        
Check if Dirichlet BCs should be applied to this node. Use the \p TOLERANCE from \p mesh_common.h as tolerance. Here a pressure value is applied if the z-coord. is equal to 4, which corresponds to one end of the pipe-mesh in this directory.
              const Real z_coo = 4.;
        
              if (fabs(curr_node(2)-z_coo) < TOLERANCE)
                {
The global number of the respective degree of freedom.
                  unsigned int dn = curr_node.dof_number(0,0,0);
        
The penalty parameter.
                  const Real penalty = 1.e10;
        
Here we apply sinusoidal pressure values for 0
                  Real p_value;
                  if (t_system.time < .002 )
                    p_value = sin(2*pi*t_system.time/.002);
                  else
                    p_value = .0;
        
Now add the contributions to the matrix and the rhs.
                  rhs.add(dn, p_value*penalty);
        
Add the panalty parameter to the global matrix if desired.
                  if (do_for_matrix)
                    matrix.add(dn, dn, penalty);
                }
            } // loop n_cnt
        }



The source file transient_ex2.C without comments:

 
  
  #include <iostream>
  #include <fstream>
  #include <algorithm>
  #include <stdio.h>
  #include <math.h>
  
  #include "libmesh/libmesh.h"
  #include "libmesh/serial_mesh.h"
  #include "libmesh/gmv_io.h"
  #include "libmesh/vtk_io.h"
  #include "libmesh/newmark_system.h"
  #include "libmesh/equation_systems.h"
  
  #include "libmesh/fe.h"
  
  #include "libmesh/quadrature_gauss.h"
  
  #include "libmesh/dense_matrix.h"
  #include "libmesh/dense_vector.h"
  
  #include "libmesh/sparse_matrix.h"
  #include "libmesh/numeric_vector.h"
  
  #include "libmesh/dof_map.h"
  
  #include "libmesh/node.h"
  
  #include "libmesh/elem.h"
  
  #include "libmesh/mesh_data.h"
  
  using namespace libMesh;
  
  void assemble_wave(EquationSystems& es,
                     const std::string& system_name);
  
  
  void apply_initial(EquationSystems& es,
                     const std::string& system_name);
  
  void fill_dirichlet_bc(EquationSystems& es,
                         const std::string& system_name);
  
  int main (int argc, char** argv)
  {
    LibMeshInit init (argc, argv);
  
    if (argc < 2)
      {
        if (libMesh::processor_id() == 0)
          std::cerr << "Usage: " << argv[0] << " [meshfile]"
                    << std::endl;
  
        libmesh_error();
      }
  
    else
      {
        std::cout << "Running " << argv[0];
  
        for (int i=1; i<argc; i++)
          std::cout << " " << argv[i];
  
        std::cout << std::endl << std::endl;
  
      }
  
    if (libMesh::default_solver_package() == LASPACK_SOLVERS)
      {
        std::cout << "WARNING! It appears you are using the\n"
                  << "LasPack solvers.  This example may not converge\n"
                  << "using LasPack, but should work OK with PETSc.\n"
                  << "http://www.mcs.anl.gov/petsc/\n"
                  << std::endl;
      }
    else if (libMesh::default_solver_package() == TRILINOS_SOLVERS)
      {
        std::cout << "WARNING! It appears you are using the\n"
                  << "Trilinos solvers.  The current libMesh Epetra\n"
                  << "interface does not allow sparse matrix addition,\n"
                  << "as is needed in this problem.  We recommend\n"
                  << "using PETSc: http://www.mcs.anl.gov/petsc/\n"
                  << std::endl;
        return 0;
      }
  
    std::string mesh_file = argv[1];
    std::cout << "Mesh file is: " << mesh_file << std::endl;
  
    libmesh_example_assert(3 <= LIBMESH_DIM, "3D support");
  
    SerialMesh mesh(init.comm());
    MeshData mesh_data(mesh);
  
    mesh.read(mesh_file, &mesh_data);
  
  
    mesh.print_info();
  
    const unsigned int result_node = 274;
  
  
    const Real delta_t = .0000625;
  
    unsigned int n_time_steps = 300;
  
    EquationSystems equation_systems (mesh);
  
    equation_systems.add_system<NewmarkSystem> ("Wave");
  
    NewmarkSystem & t_system = equation_systems.get_system<NewmarkSystem> ("Wave");
  
    t_system.add_variable("p", FIRST);
  
    t_system.attach_assemble_function  (assemble_wave);
    t_system.attach_init_function      (apply_initial);
  
    t_system.set_newmark_parameters(delta_t);
  
    equation_systems.parameters.set<Real>("speed")          = 1000.;
    equation_systems.parameters.set<Real>("fluid density")  = 1000.;
  
    t_system.time = 0.;
  
    equation_systems.init();
  
    equation_systems.print_info();
  
    std::ofstream res_out("pressure_node.res");
  
    const unsigned int res_node_no = result_node;
    const Node& res_node = mesh.node(res_node_no-1);
    unsigned int dof_no = res_node.dof_number(0,0,0);
  
    t_system.assemble();
  
    res_out << "# pressure at node " << res_node_no << "\n"
            << "# time\tpressure\n"
            << t_system.time << "\t" << 0 << std::endl;
  
  
    for (unsigned int time_step=0; time_step<n_time_steps; time_step++)
      {
        t_system.time += delta_t;
  
        t_system.update_rhs();
  
        if (time_step == 0)
          {
            equation_systems.parameters.set<bool>("Newmark set BC for Matrix") = true;
  
            fill_dirichlet_bc(equation_systems, "Wave");
  
            equation_systems.parameters.set<bool>("Newmark set BC for Matrix") = false;
          }
        else
          fill_dirichlet_bc(equation_systems, "Wave");
  
        t_system.solve();
  
        if (time_step == 30 || time_step == 60 ||
            time_step == 90 || time_step == 120 )
          {
            char buf[14];
  
  		  if (!libMesh::on_command_line("--vtk")){
  			  sprintf (buf, "out.%03d.gmv", time_step);
  
  	          GMVIO(mesh).write_equation_systems (buf,equation_systems);
  
  		  }else{
  #ifdef LIBMESH_HAVE_VTK
  			  sprintf (buf, "out_%03d.exd", time_step);
  
  	          VTKIO(mesh).write_equation_systems (buf,equation_systems);
  #endif // #ifdef LIBMESH_HAVE_VTK
  		  }
          }
  
        t_system.update_u_v_a();
  
        NumericVector<Number> &displacement
          = t_system.get_vector("displacement");
        std::vector<Number> global_displacement(displacement.size());
        displacement.localize(global_displacement);
  
        res_out << t_system.time << "\t"
                << global_displacement[dof_no]
                << std::endl;
      }
  
    return 0;
  }
  
  void assemble_wave(EquationSystems& es,
                     const std::string& system_name)
  {
    libmesh_assert_equal_to (system_name, "Wave");
  
    const MeshBase& mesh = es.get_mesh();
  
    const unsigned int dim = mesh.mesh_dimension();
  
    const Real speed = es.parameters.get<Real>("speed");
    const Real rho   = es.parameters.get<Real>("fluid density");
  
    NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);
  
    FEType fe_type = t_system.get_dof_map().variable_type(0);
  
    SparseMatrix<Number>&   stiffness = t_system.get_matrix("stiffness");
    SparseMatrix<Number>&   damping   = t_system.get_matrix("damping");
    SparseMatrix<Number>&   mass      = t_system.get_matrix("mass");
  
    NumericVector<Number>&  force     = t_system.get_vector("force");
  
    SparseMatrix<Number>&  matrix     = *t_system.matrix;
    DenseMatrix<Number>    zero_matrix;
  
    AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
  
    QGauss qrule (dim, SECOND);
  
    fe->attach_quadrature_rule (&qrule);
  
    const std::vector<Real>& JxW = fe->get_JxW();
  
    const std::vector<std::vector<Real> >& phi = fe->get_phi();
  
    const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
  
    const DofMap& dof_map = t_system.get_dof_map();
  
    DenseMatrix<Number>   Ke, Ce, Me;
    DenseVector<Number>   Fe;
  
    std::vector<dof_id_type> dof_indices;
  
    MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
    const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
  
    for ( ; el != end_el; ++el)
      {
        const Elem* elem = *el;
  
        dof_map.dof_indices (elem, dof_indices);
  
        fe->reinit (elem);
  
        {
          const unsigned int n_dof_indices = dof_indices.size();
  
          Ke.resize          (n_dof_indices, n_dof_indices);
          Ce.resize          (n_dof_indices, n_dof_indices);
          Me.resize          (n_dof_indices, n_dof_indices);
          zero_matrix.resize (n_dof_indices, n_dof_indices);
          Fe.resize          (n_dof_indices);
        }
  
        for (unsigned int qp=0; qp<qrule.n_points(); qp++)
          {
            for (unsigned int i=0; i<phi.size(); i++)
              for (unsigned int j=0; j<phi.size(); j++)
                {
                  Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
                  Me(i,j) += JxW[qp]*phi[i][qp]*phi[j][qp]
                             *1./(speed*speed);
                } // end of the matrix summation loop
          } // end of quadrature point loop
  
        {
          for (unsigned int side=0; side<elem->n_sides(); side++)
            if (!true)
              {
                AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
  
                QGauss qface(dim-1, SECOND);
  
                fe_face->attach_quadrature_rule (&qface);
  
                const std::vector<std::vector<Real> >&  phi_face = fe_face->get_phi();
  
                const std::vector<Real>& JxW_face = fe_face->get_JxW();
  
                fe_face->reinit(elem, side);
  
                const Real acc_n_value = 1.0;
  
                for (unsigned int qp=0; qp<qface.n_points(); qp++)
                  {
                    for (unsigned int i=0; i<phi_face.size(); i++)
                      {
                        Fe(i) += acc_n_value*rho
                          *phi_face[i][qp]*JxW_face[qp];
                      }
                  } // end face quadrature point loop
              } // end if (elem->neighbor(side) == NULL)
  
  
        } // end boundary condition section
  
  
        stiffness.add_matrix (Ke, dof_indices);
        damping.add_matrix   (Ce, dof_indices);
        mass.add_matrix      (Me, dof_indices);
  
        force.add_vector     (Fe, dof_indices);
  
        matrix.add_matrix(zero_matrix, dof_indices);
  
      } // end of element loop
  
    return;
  }
  
  void apply_initial(EquationSystems& es,
                     const std::string& system_name)
  {
    NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);
  
    NumericVector<Number>&  pres_vec       = t_system.get_vector("displacement");
    NumericVector<Number>&  vel_vec        = t_system.get_vector("velocity");
    NumericVector<Number>&  acc_vec        = t_system.get_vector("acceleration");
  
    pres_vec.zero();
    vel_vec.zero();
    acc_vec.zero();
  }
  
  void fill_dirichlet_bc(EquationSystems& es,
                         const std::string& system_name)
  {
    libmesh_assert_equal_to (system_name, "Wave");
  
    NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);
  
    SparseMatrix<Number>&  matrix = *t_system.matrix;
    NumericVector<Number>& rhs    = *t_system.rhs;
  
    const MeshBase& mesh = es.get_mesh();
  
    const Real pi = libMesh::pi;
  
    const bool do_for_matrix =
      es.parameters.get<bool>("Newmark set BC for Matrix");
  
    unsigned int n_nodes = mesh.n_nodes();
  
    for (unsigned int n_cnt=0; n_cnt<n_nodes; n_cnt++)
      {
        const Node& curr_node = mesh.node(n_cnt);
  
        const Real z_coo = 4.;
  
        if (fabs(curr_node(2)-z_coo) < TOLERANCE)
          {
            unsigned int dn = curr_node.dof_number(0,0,0);
  
            const Real penalty = 1.e10;
  
            Real p_value;
            if (t_system.time < .002 )
              p_value = sin(2*pi*t_system.time/.002);
            else
              p_value = .0;
  
            rhs.add(dn, p_value*penalty);
  
            if (do_for_matrix)
              matrix.add(dn, dn, penalty);
          }
      } // loop n_cnt
  }



The console output of the program:

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

Mesh file is: pipe-mesh.unv
*** Warning, This code is deprecated, and likely to be removed in future library versions! ../src/mesh/mesh_data.C, line 49, compiled Apr 19 2013 at 11:42:37 ***
 Mesh Information:
  mesh_dimension()=3
  spatial_dimension()=3
  n_nodes()=3977
    n_local_nodes()=1067
  n_elem()=3520
    n_local_elem()=880
    n_active_elem()=3520
  n_subdomains()=1
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System #0, "Wave"
    Type "Newmark"
    Variables="p" 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="FIRST", "THIRD" 
    n_dofs()=3977
    n_local_dofs()=1067
    n_constrained_dofs()=0
    n_local_constrained_dofs()=0
    n_vectors()=9
    n_matrices()=4
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 22.8167
      Average Off-Processor Bandwidth <= 1.16017
      Maximum  On-Processor Bandwidth <= 27
      Maximum Off-Processor Bandwidth <= 9
    DofMap Constraints
      Number of DoF Constraints = 0
      Number of Node Constraints = 0


 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:58:23 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=30.5408, Active time=30.2527                                                   |
 ----------------------------------------------------------------------------------------------------------------
| 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.0032      0.003154    0.0038      0.003803    0.01     0.01     |
|   build_sparsity()                 1         0.0032      0.003163    0.0111      0.011079    0.01     0.04     |
|   create_dof_constraints()         1         0.0010      0.001006    0.0010      0.001006    0.00     0.00     |
|   distribute_dofs()                1         0.0050      0.005022    0.0221      0.022149    0.02     0.07     |
|   dof_indices()                    5456      0.0253      0.000005    0.0253      0.000005    0.08     0.08     |
|   prepare_send_list()              1         0.0000      0.000027    0.0000      0.000027    0.00     0.00     |
|   reinit()                         1         0.0088      0.008766    0.0088      0.008766    0.03     0.03     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          4         0.0088      0.002195    0.0325      0.008131    0.03     0.11     |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        880       0.0054      0.000006    0.0054      0.000006    0.02     0.02     |
|   init_shape_functions()           8         0.0001      0.000015    0.0001      0.000015    0.00     0.00     |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             720       0.0023      0.000003    0.0023      0.000003    0.01     0.01     |
|   compute_map()                    160       0.0015      0.000010    0.0015      0.000010    0.01     0.01     |
|   init_reference_to_physical_map() 8         0.0002      0.000020    0.0002      0.000020    0.00     0.00     |
|                                                                                                                |
| GMVIO                                                                                                          |
|   write_nodal_data()               4         0.1061      0.026518    0.1064      0.026588    0.35     0.35     |
|                                                                                                                |
| Mesh                                                                                                           |
|   find_neighbors()                 1         0.0194      0.019424    0.0195      0.019512    0.06     0.06     |
|   read()                           1         0.0142      0.014213    0.0524      0.052437    0.05     0.17     |
|   renumber_nodes_and_elem()        2         0.0011      0.000532    0.0011      0.000532    0.00     0.00     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   broadcast()                      1         0.0066      0.006587    0.0120      0.012015    0.02     0.04     |
|   compute_hilbert_indices()        2         0.0223      0.011149    0.0223      0.011149    0.07     0.07     |
|   find_global_indices()            2         0.0027      0.001361    0.0342      0.017109    0.01     0.11     |
|   parallel_sort()                  2         0.0007      0.000339    0.0082      0.004118    0.00     0.03     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         4         0.0002      0.000042    0.1394      0.034842    0.00     0.46     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      1         0.0463      0.046274    0.0636      0.063591    0.15     0.21     |
|                                                                                                                |
| NewmarkSystem                                                                                                  |
|   initial_conditions ()            1         0.0000      0.000024    0.0000      0.000024    0.00     0.00     |
|   update_rhs ()                    300       0.1176      0.000392    0.1176      0.000392    0.39     0.39     |
|   update_u_v_a ()                  300       0.0251      0.000084    0.0251      0.000084    0.08     0.08     |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      9         0.0066      0.000729    0.0067      0.000750    0.02     0.02     |
|   broadcast()                      24        0.0013      0.000056    0.0013      0.000054    0.00     0.00     |
|   max(bool)                        4         0.0000      0.000005    0.0000      0.000005    0.00     0.00     |
|   max(scalar)                      1689      0.0068      0.000004    0.0068      0.000004    0.02     0.02     |
|   max(vector)                      342       0.0024      0.000007    0.0064      0.000019    0.01     0.02     |
|   min(bool)                        2022      0.0072      0.000004    0.0072      0.000004    0.02     0.02     |
|   min(scalar)                      1682      0.0475      0.000028    0.0475      0.000028    0.16     0.16     |
|   min(vector)                      342       0.0027      0.000008    0.0068      0.000020    0.01     0.02     |
|   probe()                          36        0.0023      0.000064    0.0023      0.000064    0.01     0.01     |
|   receive()                        36        0.0001      0.000003    0.0024      0.000067    0.00     0.01     |
|   send()                           36        0.0001      0.000002    0.0001      0.000002    0.00     0.00     |
|   send_receive()                   40        0.0002      0.000005    0.0027      0.000068    0.00     0.01     |
|   sum()                            329       0.0274      0.000083    0.0355      0.000108    0.09     0.12     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           36        0.0000      0.000001    0.0000      0.000001    0.00     0.00     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         1         0.0017      0.001701    0.0213      0.021301    0.01     0.07     |
|   set_parent_processor_ids()       1         0.0005      0.000526    0.0005      0.000526    0.00     0.00     |
|                                                                                                                |
| PetscLinearSolver                                                                                              |
|   solve()                          300       29.6619     0.098873    29.6619     0.098873    98.05    98.05    |
|                                                                                                                |
| System                                                                                                         |
|   assemble()                       1         0.0188      0.018789    0.0320      0.031996    0.06     0.11     |
|                                                                                                                |
| UNVIO                                                                                                          |
|   count_elements()                 1         0.0039      0.003902    0.0039      0.003902    0.01     0.01     |
|   count_nodes()                    1         0.0085      0.008495    0.0085      0.008495    0.03     0.03     |
|   element_in()                     1         0.0141      0.014150    0.0141      0.014150    0.05     0.05     |
|   node_in()                        1         0.0117      0.011675    0.0117      0.011675    0.04     0.04     |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            14797     30.2527                                         100.00            |
 ----------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example transient_ex2:
*  mpirun -np 4 example-devel pipe-mesh.unv -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_ex2'

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

Hosted By:
SourceForge.net Logo