/* $Id: ex8.C 2966 2008-08-10 14:28:13Z woutruijter $ *//* The Next Great Finite Element Library. *//* Copyright (C) 2003  Benjamin S. Kirk *//* This library is free software; you can redistribute it and/or *//* modify it under the terms of the GNU Lesser General Public *//* License as published by the Free Software Foundation; either *//* version 2.1 of the License, or (at your option) any later version. *//* This library is distributed in the hope that it will be useful, *//* but WITHOUT ANY WARRANTY; without even the implied warranty of *//* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU *//* Lesser General Public License for more details. *//* You should have received a copy of the GNU Lesser General Public *//* License along with this library; if not, write to the Free Software *//* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA */ // <h1>Example 8 - The Wave Equation</h1> // // 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.h"#include "mesh.h"#include "gmv_io.h"#include "vtk_io.h"#include "newmark_system.h"#include "equation_systems.h"// Define the Finite Element object.#include "fe.h"// Define Gauss quadrature rules.#include "quadrature_gauss.h"// Define useful datatypes for finite element#include "dense_matrix.h"#include "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 "sparse_matrix.h"#include "numeric_vector.h"// Define the DofMap, which handles degree of freedom// indexing.#include "dof_map.h"// The definition of a vertex associated with a Mesh.#include "node.h"// The definition of a geometric element#include "elem.h"// Defines the MeshData class, which allows you to store// data about the mesh when reading in files, etc.#include "mesh_data.h"// 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 programint main (int argc, char** argv){  // Initialize Petsc, like in example 2.  LibMeshInit init (argc, argv);#ifdef ENABLE_PARMESH  if (libMesh::processor_id() == 0)    std::cerr << "ERROR: This example directly references\n"              << "all mesh nodes and is incompatible with"              << "ParallelMesh use."              << std::endl;#else  // 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).  Print a warning to the user if PETSc  // is not available, or if they are using LasPack solvers.#ifdef HAVE_PETSC  if ((libMesh::on_command_line("--use-laspack")) ||      (libMesh::on_command_line("--disable-petsc")))#endif    {      std::cout << "WARNING! It appears you are using the\n"                << "LasPack solvers.  ex8 is known not to converge\n"                << "using LasPack, but should work OK with PETSc.\n"                << "If possible, download and install the PETSc\n"                << "library from www-unix.mcs.anl.gov/petsc/petsc-2/\n"                << std::endl;    }    // 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;  // For now, restrict to dim=3, though this  // may easily be changed, see example 4  const unsigned int dim = 3;  // Create a dim-dimensional mesh.  Mesh mesh (dim);  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.;  // Store the current time as an  // \p EquationSystems parameter, so that  // \p fill_dirichlet_bc() can access it.  equation_systems.parameters.set<Real>("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  Real t_time = 0.;  res_out << "# pressure at node " << res_node_no << "\n"          << "# time\tpressure\n"          << t_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_time += delta_t;      equation_systems.parameters.set<Real>("time")  = t_time;      // 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{			  // VTK viewers are generally not happy with two dots in a filename			  sprintf (buf, "out_%03d.pvtu", time_step);	          VTKIO(mesh).write_equation_systems (buf,equation_systems);		  }        }      // 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);