ex8.c

一个用来实现偏微分方程中网格的计算库

/* $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 program
int 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);