equation_systems_io.c

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

// $Id: equation_systems_io.C 2858 2008-06-13 18:59:07Z benkirk $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// 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


#include "libmesh_common.h"
#include "libmesh_logging.h"


// C++ Includes
#include <cstdio> // for std::sprintf

// Local Includes
#include "equation_systems.h"
#include "mesh_base.h"
#include "mesh_tools.h"
#include "parallel_mesh.h"
#include "parallel.h"
#include "serial_mesh.h"
#include "xdr_cxx.h"

// Forward Declarations

// Anonymous namespace for implementation details.
namespace {
  std::string local_file_name (const std::string &name)
  {
    std::string basename(name);
    char buf[256];

    if (basename.size() - basename.rfind(".bz2") == 4)
      {
	basename.erase(basename.end()-4, basename.end());
	std::sprintf(buf, "%s.%04d.bz2", basename.c_str(), libMesh::processor_id());
      }
    else if (basename.size() - basename.rfind(".gz") == 3)
      {
	basename.erase(basename.end()-3, basename.end());
	std::sprintf(buf, "%s.%04d.gz", basename.c_str(), libMesh::processor_id());
      }
    else
      std::sprintf(buf, "%s.%04d", basename.c_str(), libMesh::processor_id());
      
    return std::string(buf);
  }
}




// ------------------------------------------------------------
// EquationSystem class implementation
void EquationSystems::read (const std::string& name,
			    const libMeshEnums::XdrMODE mode,
                            const unsigned int read_flags)
{
  /**
   * This program implements the output of an 
   * EquationSystems object.  This warrants some 
   * documentation.  The output file essentially
   * consists of 11 sections:
     \verbatim
     1.) A version header (for non-'legacy' formats, libMesh-0.7.0 and greater). 
     2.) The number of individual equation systems (unsigned int)
     
       for each system
                                                          
        3.)  The name of the system (string)
        4.)  The type of the system (string)
    
        handled through System::read():
    
     +-------------------------------------------------------------+
     |  5.) The number of variables in the system (unsigned int)   |
     |                                                             |
     |   for each variable in the system                           |
     |                                                             |
     |    6.) The name of the variable (string)                    |
     |                                                             |
     |    7.) Combined in an FEType:                               |
     |         - The approximation order(s) of the variable (Order |
     |           Enum, cast to int/s)                              |
     |         - The finite element family/ies of the variable     |
     |           (FEFamily Enum, cast to int/s)                    |
     |                                                             |
     |   end variable loop                                         |
     |                                                             |
     | 8.) The number of additional vectors (unsigned int),        |
     |                                                             |
     |    for each additional vector in the equation system object |
     |                                                             |
     |    9.) the name of the additional vector  (string)          |
     +-------------------------------------------------------------+
    
     end system loop
    
    
       for each system, handled through System::read_{serialized,parallel}_data():
       
     +--------------------------------------------------------------+
     | 10.) The global solution vector, re-ordered to be node-major |
     |     (More on this later.)                                    |
     |                                                              |
     |    for each additional vector in the equation system object  |
     |                                                              |
     |    11.) The global additional vector, re-ordered to be       |
     |         node-major (More on this later.)                     |
     +--------------------------------------------------------------+
    
     end system loop
   \endverbatim
   *
   * Note that the actual IO is handled through the Xdr class 
   * (to be renamed later?) which provides a uniform interface to 
   * both the XDR (eXternal Data Representation) interface and standard
   * ASCII output.  Thus this one section of code will read XDR or ASCII
   * files with no changes.
   */

   // Set booleans from the read_flags argument
   const bool read_header          = read_flags & EquationSystems::READ_HEADER;
   const bool read_data            = read_flags & EquationSystems::READ_DATA;
   const bool read_additional_data = read_flags & EquationSystems::READ_ADDITIONAL_DATA;
   const bool read_legacy_format   = read_flags & EquationSystems::READ_LEGACY_FORMAT;
         bool read_parallel_files  = false;

  // This will unzip a file with .bz2 as the extension, otherwise it
  // simply returns the name if the file need not be unzipped.
  Xdr io ((libMesh::processor_id() == 0) ? name : "", mode);
  libmesh_assert (io.reading());
	  
  {
    // 1.)
    // Read the version header.
    std::string version = "legacy";
    if (!read_legacy_format)
      {
	if (libMesh::processor_id() == 0) io.data(version);	
	Parallel::broadcast(version);
	
	// All processors have the version header, if it does not contain
	// "libMesh" something then it is a legacy file.
	if (!(version.find("libMesh") < version.size()))
	  {
	    io.close();
	    
	    // Recursively call this read() function but with the 
	    // EquationSystems::READ_LEGACY_FORMAT bit set.
	    this->read (name, mode, (read_flags | EquationSystems::READ_LEGACY_FORMAT));
	    return;
	  }

	read_parallel_files = (version.rfind(" parallel") < version.size());
      }
    else
      deprecated();

  START_LOG("read()","EquationSystems");
      
    // 2.)  
    // Read the number of equation systems
    unsigned int n_sys=0;
    if (libMesh::processor_id() == 0) io.data (n_sys);
    Parallel::broadcast(n_sys);

    for (unsigned int sys=0; sys<n_sys; sys++)
      {
	// 3.)
	// Read the name of the sys-th equation system
	std::string sys_name;      
	if (libMesh::processor_id() == 0) io.data (sys_name);
	Parallel::broadcast(sys_name);
	
	// 4.)
	// Read the type of the sys-th equation system
	std::string sys_type;	
	if (libMesh::processor_id() == 0) io.data (sys_type);
	Parallel::broadcast(sys_type);
	
	if (read_header)
	  this->add_system (sys_type, sys_name);

	// 5.) - 9.)
	// Let System::read_header() do the job
	System& new_system = this->get_system(sys_name);	  
	new_system.read_header (io,
				read_header,
				read_additional_data,
				read_legacy_format);
      }
  }
      


  // Now we are ready to initialize the underlying data
  // structures. This will initialize the vectors for 
  // storage, the dof_map, etc...
  if (read_header) 
    this->init();

  // 10.) & 11.)
  // Read and set the numeric vector values
  if (read_data)
    {
      // the EquationSystems::read() method should look constant from the mesh
      // perspective, but we need to assign a temporary numbering to the nodes
      // and elements in the mesh, which requires that we abuse const_cast
      if (!read_legacy_format)
	{
	  MeshBase &mesh = const_cast<MeshBase&>(this->get_mesh());
	  MeshTools::Private::globally_renumber_nodes_and_elements(mesh);
	}

      Xdr local_io (read_parallel_files ? local_file_name(name) : "", mode);
   
      std::map<std::string, System*>::iterator
	pos = _systems.begin();