system_io.c

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

// $Id: system_io.C 2955 2008-08-01 20:44:59Z roystgnr $

// 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 "parallel.h"

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

// Local Includes
#include "system.h"
#include "mesh.h"
#include "mesh_tools.h"
#include "elem.h"
#include "xdr_cxx.h"
#include "numeric_vector.h"



// Anonymous namespace for implementation details.
namespace {
  
  // Comments:
  // ---------
  // - The io_blksize governs how many nodes or elements will be treated as a single block
  //   when performing parallel IO.
  // - This parameter only loosely affects the size of the actual IO buffer as this depends
  //   on the number of components a given variable has for the nodes/elements in the block.
  // - When reading/writing each processor uses an ID map which is 3*io_blksize*sizeof(unsigned int) bytes
  //   long, so if io_blksize=256000 we would expect that buffer alone to be ~3Mb.
  // - In general, an increase in io_blksize should increase the efficiency of the parallel
  //   read/writes by reducing the number of MPI messages at the expense of memory.
  // - If the library exhausts memory during IO you might reduce this parameter.
  
  const unsigned int io_blksize = 256000;


  /**
   * Comparison object to use with DofObject pointers.  This sorts by id(),
   * so when we iterate over a set of DofObjects we visit the objects in
   * order of increasing ID.
   */
  struct CompareDofObjectsByID
  {
    bool operator()(const DofObject *a,
		    const DofObject *b) const
    {
      libmesh_assert (a);
      libmesh_assert (b);

      return a->id() < b->id();
    }
  };
}



// ------------------------------------------------------------
// System class implementation
void System::read_header (Xdr& io,
			  const bool read_header,
			  const bool read_additional_data,
			  const bool read_legacy_format)
{
  // This method implements the input of a
  // System object, embedded in the output of
  // an EquationSystems<T_sys>.  This warrants some 
  // documentation.  The output file essentially
  // consists of 5 sections:
  //
  // for this system
  //
  //   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 system object
  // 
  //     9.) the name of the additional vector  (string)
  //
  // end system
  libmesh_assert (io.reading());
  
  // Possibly clear data structures and start from scratch.
  if (read_header)
    this->clear ();
  
  {
    // 5.) 
    // Read the number of variables in the system
    unsigned int n_vars=0;
    if (libMesh::processor_id() == 0) io.data (n_vars);
    Parallel::broadcast(n_vars);
      
    for (unsigned int var=0; var<n_vars; var++)
      {	            
	// 6.)
	// Read the name of the var-th variable
	std::string var_name;	  
	if (libMesh::processor_id() == 0) io.data (var_name);
	Parallel::broadcast(var_name);
	      	
	// 7.)
	// Read the approximation order(s) of the var-th variable 
	int order=0;	  
	if (libMesh::processor_id() == 0) io.data (order);
	Parallel::broadcast(order);
	
#ifdef ENABLE_INFINITE_ELEMENTS
	
	// do the same for radial_order
	int rad_order=0;	  
	if (libMesh::processor_id() == 0) io.data(rad_order);
	Parallel::broadcast(rad_order);

#endif

	// Read the finite element type of the var-th variable 
	int fam=0;
	if (libMesh::processor_id() == 0) io.data (fam);
	Parallel::broadcast(fam);
	FEType type;
	type.order  = static_cast<Order>(order);
	type.family = static_cast<FEFamily>(fam);

	// Check for incompatibilities.  The shape function indexing was
	// changed for the monomial and xyz finite element families to 
	// simplify extension to arbitrary p.  The consequence is that 
	// old restart files will not be read correctly.  This is expected
	// to be an unlikely occurance, but catch it anyway.
	if (read_legacy_format)
	  if ((type.family == MONOMIAL || type.family == XYZ) && 
	      ((type.order > 2 && this->get_mesh().mesh_dimension() == 2) ||
	       (type.order > 1 && this->get_mesh().mesh_dimension() == 3)))
	    {
	      here();
	      std::cout << "*****************************************************************\n"
			<< "* WARNING: reading a potentially incompatible restart file!!!   *\n"
			<< "*  contact [e-mail] for more details *\n"
			<< "*****************************************************************"
			<< std::endl;
	    }
			
#ifdef ENABLE_INFINITE_ELEMENTS
	
	// Read additional information for infinite elements	
	int radial_fam=0;
	int i_map=0;
	
	if (libMesh::processor_id() == 0) io.data (radial_fam);
	Parallel::broadcast(radial_fam);
	if (libMesh::processor_id() == 0) io.data (i_map);
	Parallel::broadcast(i_map);
	
	type.radial_order  = static_cast<Order>(rad_order);
	type.radial_family = static_cast<FEFamily>(radial_fam);
	type.inf_map       = static_cast<InfMapType>(i_map);	  

#endif

	if (read_header) 
	  this->add_variable (var_name, type);
      }
  }

  // 8.)  
  // Read the number of additional vectors.  
  unsigned int n_vectors=0;  
  if (libMesh::processor_id() == 0) io.data (n_vectors);
  Parallel::broadcast(n_vectors);
  
  // If n_vectors > 0, this means that write_additional_data
  // was true when this file was written.  We will need to
  // make use of this fact later.
  if (n_vectors > 0)
    this->_additional_data_written = true;  
  
  for (unsigned int vec=0; vec<n_vectors; vec++)
    {
      // 9.)
      // Read the name of the vec-th additional vector
      std::string vec_name;      
      if (libMesh::processor_id() == 0) io.data (vec_name);
      Parallel::broadcast(vec_name);
      
      if (read_additional_data)
	{
	  // sanity checks
	  libmesh_assert(this->_can_add_vectors);
	  // Some systems may have added their own vectors already
//	  libmesh_assert(this->_vectors.count(vec_name) == 0);

	  this->add_vector(vec_name);
	}
    }
}



void System::read_legacy_data (Xdr& io,
			       const bool read_additional_data)
{
  deprecated();
  
  // This method implements the output of the vectors
  // contained in this System object, embedded in the 
  // output of an EquationSystems<T_sys>. 
  //
  //   10.) The global solution vector, re-ordered to be node-major 
  //       (More on this later.)                                    
  //                                                                
  //      for each additional vector in the object          
  //                                                                
  //      11.) The global additional vector, re-ordered to be       
  //           node-major (More on this later.)
  libmesh_assert (io.reading());

  // read and reordering buffers
  std::vector<Number> global_vector;
  std::vector<Number> reordered_vector;

  // 10.)
  // Read and set the solution vector
  {	
    if (libMesh::processor_id() == 0) io.data (global_vector);	  
    Parallel::broadcast(global_vector);
    
    // Remember that the stored vector is node-major.
    // We need to put it into whatever application-specific
    // ordering we may have using the dof_map.
    reordered_vector.resize(global_vector.size());

    //std::cout << "global_vector.size()=" << global_vector.size() << std::endl;
    //std::cout << "this->n_dofs()=" << this->n_dofs() << std::endl;
    
    libmesh_assert (global_vector.size() == this->n_dofs());
	
    unsigned int cnt=0;

    const unsigned int sys     = this->number();
    const unsigned int n_vars  = this->n_vars();

    for (unsigned int var=0; var<n_vars; var++)
      {	
	// First reorder the nodal DOF values
	{
	  MeshBase::node_iterator
	    it  = this->get_mesh().nodes_begin(),
	    end = this->get_mesh().nodes_end();
	
  	  for (; it != end; ++it)
	    for (unsigned int index=0; index<(*it)->n_comp(sys,var); index++)
	      {
	        libmesh_assert ((*it)->dof_number(sys, var, index) !=
			DofObject::invalid_id);
		
                libmesh_assert (cnt < global_vector.size());
                
                reordered_vector[(*it)->dof_number(sys, var, index)] =
	  	  global_vector[cnt++]; 
	      }
	}
	
	// Then reorder the element DOF values
	{
	  MeshBase::element_iterator
	    it  = this->get_mesh().active_elements_begin(),
	    end = this->get_mesh().active_elements_end();

	  for (; it != end; ++it)
	    for (unsigned int index=0; index<(*it)->n_comp(sys,var); index++)
	      {  
		libmesh_assert ((*it)->dof_number(sys, var, index) !=
			DofObject::invalid_id);
		
		libmesh_assert (cnt < global_vector.size());
		
		reordered_vector[(*it)->dof_number(sys, var, index)] =
		  global_vector[cnt++]; 
	      }
	}
      }
	    
    *(this->solution) = reordered_vector;
  }

  // For each additional vector, simply go through the list.
  // ONLY attempt to do this IF additional data was actually
  // written to the file for this system (controlled by the
  // _additional_data_written flag).  
  if (this->_additional_data_written)
    {
      std::map<std::string, NumericVector<Number>* >::iterator
	pos = this->_vectors.begin();
  
      for (; pos != this->_vectors.end(); ++pos)
	{
	  // 11.) 	   
	  // Read the values of the vec-th additional vector.
	  // Prior do _not_ clear, but fill with zero, since the
	  // additional vectors _have_ to have the same size
	  // as the solution vector
	  std::fill (global_vector.begin(), global_vector.end(), libMesh::zero);

	  if (libMesh::processor_id() == 0) io.data (global_vector);
	  Parallel::broadcast(global_vector);

	  // If read_additional_data==true, then we will keep this vector, otherwise
	  // we are going to throw it away.
	  if (read_additional_data)
	    {
	      // Remember that the stored vector is node-major.
	      // We need to put it into whatever application-specific
	      // ordering we may have using the dof_map.
	      std::fill (reordered_vector.begin(),
			 reordered_vector.end(),
			 libMesh::zero);
	
	      reordered_vector.resize(global_vector.size());	

	      libmesh_assert (global_vector.size() == this->n_dofs());
	
	      unsigned int cnt=0;

	      const unsigned int sys     = this->number();
	      const unsigned int n_vars  = this->n_vars();
	
	      for (unsigned int var=0; var<n_vars; var++)
		{
		  // First reorder the nodal DOF values
		  {
		    MeshBase::node_iterator
		      it  = this->get_mesh().nodes_begin(),
		      end = this->get_mesh().nodes_end();

		    for (; it!=end; ++it)
		      for (unsigned int index=0; index<(*it)->n_comp(sys,var); index++)
			{
			  libmesh_assert ((*it)->dof_number(sys, var, index) !=
				  DofObject::invalid_id);
			  
			  libmesh_assert (cnt < global_vector.size());
			  
			  reordered_vector[(*it)->dof_number(sys, var, index)] =
			    global_vector[cnt++]; 
			}
		  }

		  // Then reorder the element DOF values
		  {
		    MeshBase::element_iterator
		      it  = this->get_mesh().active_elements_begin(),
		      end = this->get_mesh().active_elements_end();

		    for (; it!=end; ++it)
		      for (unsigned int index=0; index<(*it)->n_comp(sys,var); index++)
			{  
			  libmesh_assert ((*it)->dof_number(sys, var, index) !=
				  DofObject::invalid_id);
			  
			  libmesh_assert (cnt < global_vector.size());
			  
			  reordered_vector[(*it)->dof_number(sys, var, index)] =
			    global_vector[cnt++]; 
			}
		  }
		}
	    
	      // use the overloaded operator=(std::vector) to assign the values
	      *(pos->second) = reordered_vector;
	    }
	}
    } // end if (_additional_data_written)    
}