ex6.c

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

/* $Id: ex6.C 2837 2008-05-08 17:23:37Z roystgnr $ */

/* 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 6 - Infinite Elements for the Wave Equation</h1>
 //
 // This is the sixth example program.  It builds on
 // the previous examples, and introduces the Infinite
 // Element class.  Note that the library must be compiled
 // with Infinite Elements enabled.  Otherwise, this
 // example will abort.
 // This example intends to demonstrate the similarities
 // between the \p FE and the \p InfFE classes in libMesh.
 // The matrices are assembled according to the wave equation.
 // However, for practical applications a time integration
 // scheme (as introduced in subsequent examples) should be
 // used.

// C++ include files that we need
#include <iostream>
#include <algorithm>
#include <math.h>

// Basic include file needed for the mesh functionality.
#include "libmesh.h"
#include "mesh.h"
#include "mesh_generation.h"
#include "gmv_io.h"
#include "linear_implicit_system.h"
#include "equation_systems.h"

// Define the Finite and Infinite Element object.
#include "fe.h"
#include "inf_fe.h"
#include "inf_elem_builder.h"

// Define Gauss quadrature rules.
#include "quadrature_gauss.h"

// Define useful datatypes for finite element
// matrix and vector components.
#include "sparse_matrix.h"
#include "numeric_vector.h"
#include "dense_matrix.h"
#include "dense_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"

// Function prototype.  This is similar to the Poisson
// assemble function of example 4.  
void assemble_wave (EquationSystems& es,
                    const std::string& system_name);

// Begin the main program.
int main (int argc, char** argv)
{
  // Initialize libMesh, like in example 2.
  LibMeshInit init (argc, argv);
  
  // This example requires Infinite Elements   
#ifndef ENABLE_INFINITE_ELEMENTS
  if (libMesh::processor_id() == 0)
    std::cerr << "ERROR: This example requires the library to be" << std::endl
              << "compiled with Infinite Element support!" << std::endl;

  return 0;

#else
  
  // For the moment, only allow 3D     
  const unsigned int dim = 3; 
  
  // Tell the user what we are doing.
  std::cout << "Running ex6 with dim = " << dim << std::endl << std::endl;        
  
  // Create a mesh with user-defined dimension 
  Mesh mesh (dim);

  // Use the internal mesh generator to create elements
  // on the square [-1,1]^3, of type Hex8.
  MeshTools::Generation::build_cube (mesh,
                                     4, 4, 4,
                                     -1., 1.,
                                     -1., 1.,
                                     -1., 1.,
                                     HEX8);
  
  // Print information about the mesh to the screen.
  mesh.print_info();

  // Write the mesh before the infinite elements are added
  GMVIO(mesh).write ("orig_mesh.gmv");

  // Normally, when a mesh is imported or created in
  // libMesh, only conventional elements exist.  The infinite
  // elements used here, however, require prescribed
  // nodal locations (with specified distances from an imaginary
  // origin) and configurations that a conventional mesh creator 
  // in general does not offer.  Therefore, an efficient method
  // for building infinite elements is offered.  It can account
  // for symmetry planes and creates infinite elements in a fully
  // automatic way.
  //
  // Right now, the simplified interface is used, automatically
  // determining the origin.  Check \p MeshBase for a generalized
  // method that can even return the element faces of interior
  // vibrating surfaces.  The \p bool determines whether to be 
  // verbose.
  InfElemBuilder builder(mesh);
  builder.build_inf_elem(true);

  // Print information about the mesh to the screen.
  mesh.print_info();

  // Write the mesh with the infinite elements added.
  // Compare this to the original mesh.
  GMVIO(mesh).write ("ifems_added.gmv");

  // After building infinite elements, we have to let 
  // the elements find their neighbors again.
  mesh.find_neighbors();
  
  // Create an equation systems object, where \p ThinSystem
  // offers only the crucial functionality for solving a 
  // system.  Use \p ThinSystem when you want the sleekest
  // system possible.
  EquationSystems equation_systems (mesh);
  
  // Declare the system and its variables.
  // Create a system named "Wave".  This can
  // be a simple, steady system
  equation_systems.add_system<LinearImplicitSystem> ("Wave");
        
  // Create an FEType describing the approximation
  // characteristics of the InfFE object.  Note that
  // the constructor automatically defaults to some
  // sensible values.  But use \p FIRST order 
  // approximation.
  FEType fe_type(FIRST);
  
  // Add the variable "p" to "Wave".  Note that there exist
  // various approaches in adding variables.  In example 3, 
  // \p add_variable took the order of approximation and used
  // default values for the \p FEFamily, while here the \p FEType 
  // is used.
  equation_systems.get_system("Wave").add_variable("p", fe_type);
  
  // Give the system a pointer to the matrix assembly
  // function.
  equation_systems.get_system("Wave").attach_assemble_function (assemble_wave);
  
  // 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")          = 1.;
  equation_systems.parameters.set<Real>("fluid density")  = 1.;
  
  // Initialize the data structures for the equation system.
  equation_systems.init();
  
  // Prints information about the system to the screen.
  equation_systems.print_info();

  // Solve the system "Wave".
  equation_systems.get_system("Wave").solve();
  
  // Write the whole EquationSystems object to file.
  // For infinite elements, the concept of nodal_soln()
  // is not applicable. Therefore, writing the mesh in
  // some format @e always gives all-zero results at
  // the nodes of the infinite elements.  Instead,
  // use the FEInterface::compute_data() methods to
  // determine physically correct results within an
  // infinite element.
  equation_systems.write ("eqn_sys.dat", libMeshEnums::WRITE);
  
  // All done.  
  return 0;

#endif // else part of ifndef ENABLE_INFINITE_ELEMENTS
}

// This function assembles the system matrix and right-hand-side
// for the discrete form of 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 (system_name == "Wave");


#ifdef ENABLE_INFINITE_ELEMENTS
  
  // Get a constant reference to the mesh object.
  const MeshBase& mesh = es.get_mesh();

  // Get a reference to the system we are solving.
  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Wave");
  
  // 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 = system.get_dof_map();
  
  // 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");
  
  // Get a constant reference to the Finite Element type
  // for the first (and only) variable in the system.
  const FEType& fe_type = dof_map.variable_type(0);
  
  // Build a Finite Element object of the specified type.  Since the
  // \p FEBase::build() member dynamically creates memory we will
  // store the object as an \p AutoPtr<FEBase>.  Check ex5 for details.
  AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
  
  // Do the same for an infinite element.
  AutoPtr<FEBase> inf_fe (FEBase::build_InfFE(dim, fe_type));