<commit_before>// The libMesh Finite Element Library.
// Copyright (C) 2002-2020 <NAME>, <NAME>, <NAME>

// 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>Transient Example 2 - The Newmark System and the Wave Equation</h1>
// \author <NAME>
// \date 2003
//
// 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 <sstream>
#include <algorithm>
#include <stdio.h>
#include <math.h>

// Basic include file needed for the mesh functionality.
#include "libmesh/libmesh.h"
#include "libmesh/replicated_mesh.h"
#include "libmesh/gmv_io.h"
#include "libmesh/vtk_io.h"
#include "libmesh/newmark_system.h"
#include "libmesh/equation_systems.h"
#include "libmesh/enum_solver_package.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 and Eigen 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"

// 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);

  // This example requires a linear solver package.
  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
                           "--enable-petsc, --enable-trilinos, or --enable-eigen");

  // Check for proper usage.
  libmesh_error_msg_if(argc < 2, "Usage: " << argv[0] << " [meshfile]");

  // Tell the user what we are doing.
  libMesh::out << "Running " << argv[0];

  for (int i=1; i<argc; i++)
    libMesh::out << " " << argv[i];

  libMesh::out << std::endl << std::endl;

  // LasPack solvers don't work so well for this example, Trilinos doesn't work at all.
  // PETSc and Eigen both work...
  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS || \
                           libMesh::default_solver_package() == EIGEN_SOLVERS, "--enable-petsc");

  // Get the name of the mesh file
  // from the command line.
  std::string mesh_file = argv[1];
  libMesh::out << "Mesh file is: " << mesh_file << std::endl;

  // Skip this 3D example if libMesh was compiled as 1D or 2D-only.
  libmesh_example_requires(3 <= LIBMESH_DIM, "3D support");

  // Create a mesh.
  // This example directly references all mesh nodes and is
  // incompatible with DistributedMesh use.
  //
  // Create a ReplicatedMesh object, with dimension to be overridden
  // later, distributed across the default MPI communicator.
  ReplicatedMesh mesh(init.comm());

  // Read the meshfile specified on the command line.
  mesh.read(mesh_file);

  // 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 NewmarkSystem can
  // compute integration constants.  Here we simply use
  // pass only the time step and use default values
  // for alpha=.25  and delta=.5.
  t_system.set_newmark_parameters(delta_t);

  // Set the speed of sound and fluid density
  // as EquationSystems parameter,
  // so that 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_ref(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 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
      // 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 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)
        {
          std::ostringstream file_name;

#ifdef LIBMESH_HAVE_VTK
          file_name << "out_"
                    << std::setw(3)
                    << std::setfill('0')
                    << std::right
                    << time_step
                    << ".pvtu";

          VTKIO(mesh).write_equation_systems (file_name.str(), equation_systems);
#else

          file_name << "out."
                    << std::setw(3)
                    << std::setfill('0')
                    << std::right
                    << time_step
                    << ".gmv";

          GMVIO(mesh).write_equation_systems (file_name.str(), equation_systems);
#endif
        }

      // 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;<commit_msg>Comment out now-clearly-unused variable too<commit_after>// The libMesh Finite Element Library.
// Copyright (C) 2002-2020 <NAME>, <NAME>, <NAME>

// 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>Transient Example 2 - The Newmark System and the Wave Equation</h1>
// \author <NAME>
// \date 2003
//
// 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 <sstream>
#include <algorithm>
#include <stdio.h>
#include <math.h>

// Basic include file needed for the mesh functionality.
#include "libmesh/libmesh.h"
#include "libmesh/replicated_mesh.h"
#include "libmesh/gmv_io.h"
#include "libmesh/vtk_io.h"
#include "libmesh/newmark_system.h"
#include "libmesh/equation_systems.h"
#include "libmesh/enum_solver_package.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 and Eigen 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"

// 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);

  // This example requires a linear solver package.
  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
                           "--enable-petsc, --enable-trilinos, or --enable-eigen");

  // Check for proper usage.
  libmesh_error_msg_if(argc < 2, "Usage: " << argv[0] << " [meshfile]");

  // Tell the user what we are doing.
  libMesh::out << "Running " << argv[0];

  for (int i=1; i<argc; i++)
    libMesh::out << " " << argv[i];

  libMesh::out << std::endl << std::endl;

  // LasPack solvers don't work so well for this example, Trilinos doesn't work at all.
  // PETSc and Eigen both work...
  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS || \
                           libMesh::default_solver_package() == EIGEN_SOLVERS, "--enable-petsc");

  // Get the name of the mesh file
  // from the command line.
  std::string mesh_file = argv[1];
  libMesh::out << "Mesh file is: " << mesh_file << std::endl;

  // Skip this 3D example if libMesh was compiled as 1D or 2D-only.
  libmesh_example_requires(3 <= LIBMESH_DIM, "3D support");

  // Create a mesh.
  // This example directly references all mesh nodes and is
  // incompatible with DistributedMesh use.
  //
  // Create a ReplicatedMesh object, with dimension to be overridden
  // later, distributed across the default MPI communicator.
  ReplicatedMesh mesh(init.comm());

  // Read the meshfile specified on the command line.
  mesh.read(mesh_file);

  // 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 NewmarkSystem can
  // compute integration constants.  Here we simply use
  // pass only the time step and use default values
  // for alpha=.25  and delta=.5.
  t_system.set_newmark_parameters(delta_t);

  // Set the speed of sound and fluid density
  // as EquationSystems parameter,
  // so that 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_ref(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 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
      // 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 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)
        {
          std::ostringstream file_name;

#ifdef LIBMESH_HAVE_VTK
          file_name << "out_"
                    << std::setw(3)
                    << std::setfill('0')
                    << std::right
                    << time_step
                    << ".pvtu";

          VTKIO(mesh).write_equation_systems (file_name.str(), equation_systems);
#else

          file_name << "out."
                    << std::setw(3)
                    << std::setfill('0')
                    << std::right
                    << time_step
                    << ".gmv";

          GMVIO(mesh).write_equation_systems (file_name.str(), equation_systems);
#endif
        }

      // 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 to a local variable.
  const Real speed = es.parameters.get<Real>("speed");

  // If we added Neumann conditions we would need density too
  // 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;
{'fraction_non_alphanumeric': 0.10031790886612504, 'fraction_numerical': 0.007457121551081283, 'mean_word_length': 2.8617762958472266, 'pattern_counts': {'":': 0, '<': 207, '<?xml version=': 0, '>': 59, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5706280', 'n_tokens_mistral': 7657, 'n_tokens_neox': 7164, 'n_words': 2708}