xdr_io.c

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

// $Id: legacy_xdr_io.C 2560 2007-12-03 17:52:20Z 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



// C++ includes
#include <iostream>
#include <iomanip>

#include <vector>
#include <string>

// Local includes
#include "xdr_io.h"
#include "legacy_xdr_io.h"
#include "xdr_cxx.h"
#include "enum_xdr_mode.h"
#include "mesh_base.h"
#include "node.h"
#include "elem.h"
#include "boundary_info.h"
#include "parallel.h"
#include "mesh_tools.h"
#include "partitioner.h"
#include "libmesh_logging.h"


//-----------------------------------------------
// anonymous namespace for implementation details
namespace {
  struct ElemBCData
  {
    unsigned int       elem_id;
    unsigned short int side;
    short int          bc_id;
    
    // Default constructor
    ElemBCData (unsigned int       elem_id_in=0,
		unsigned short int side_in=0,
		short int          bc_id_in=0) :
      elem_id(elem_id_in),
      side(side_in),
      bc_id(bc_id_in)
    {}

    // comparison operator
    bool operator < (const ElemBCData &other) const
    {
      if (this->elem_id == other.elem_id)
	return (this->side < other.side);
      
      return this->elem_id < other.elem_id;
    }  
  };

  // comparison operator
  bool operator < (const unsigned int &other_elem_id,
		   const ElemBCData &elem_bc)
  {
    return other_elem_id < elem_bc.elem_id;
  }

  bool operator < (const ElemBCData &elem_bc,
		   const unsigned int &other_elem_id)
    
  {
    return elem_bc.elem_id < other_elem_id;
  }


  // For some reason SunStudio does not seem to accept the above
  // comparison functions for use in
  // std::equal_range (ElemBCData::iterator, ElemBCData::iterator, unsigned int);
  struct CompareIntElemBCData
  {
    bool operator()(const unsigned int &other_elem_id,
		    const ElemBCData &elem_bc)
    {
      return other_elem_id < elem_bc.elem_id;
    }

    bool operator()(const ElemBCData &elem_bc,
		    const unsigned int &other_elem_id)    
    {
      return elem_bc.elem_id < other_elem_id;
    }
  };
}



// ------------------------------------------------------------
// XdrIO static data
const unsigned int XdrIO::io_blksize = 128000;



// ------------------------------------------------------------
// XdrIO members
XdrIO::XdrIO (MeshBase& mesh, const bool binary) :
  MeshInput<MeshBase> (mesh,/* is_parallel_format = */ true),
  MeshOutput<MeshBase>(mesh,/* is_parallel_format = */ true),
  _binary             (binary),
  _legacy             (false),
  _version            ("libMesh-0.7.0+"),
  _bc_file_name       ("n/a"),
  _partition_map_file ("n/a"),
  _subdomain_map_file ("n/a"),
  _p_level_file       ("n/a")
{
}



XdrIO::XdrIO (const MeshBase& mesh, const bool binary) :
  MeshOutput<MeshBase>(mesh,/* is_parallel_format = */ true),
  _binary (binary)
{
}



XdrIO::~XdrIO ()
{
}



void XdrIO::write (const std::string& name)
{
  if (this->legacy())
    {
      deprecated();
      LegacyXdrIO(MeshOutput<MeshBase>::mesh(), this->binary()).write(name);
      return;
    }
  
  Xdr io ((libMesh::processor_id() == 0) ? name : "", this->binary() ? ENCODE : WRITE);

  START_LOG("write()","XdrIO");
  
  // convenient reference to our mesh
  const MeshBase &mesh = MeshOutput<MeshBase>::mesh();
  
  unsigned int
    n_elem     = mesh.n_elem(),
    n_nodes    = mesh.n_nodes(),
    n_bcs      = mesh.boundary_info->n_boundary_conds(),
    n_p_levels = MeshTools::n_p_levels (mesh); 


  //-------------------------------------------------------------
  // For all the optional files -- the default file name is "n/a".
  // However, the user may specify an optional external file. 
  
  // If there are BCs and the user has not already provided a
  // file name then write to "."
  if (n_bcs &&
      this->boundary_condition_file_name() == "n/a")
    this->boundary_condition_file_name() = ".";
  
  // If there are more than one subdomains and the user has not specified an 
  // external file then write the subdomain mapping to the default file "."  
  if ((mesh.n_subdomains() > 1) && 
      (this->subdomain_map_file_name() == "n/a"))
    this->subdomain_map_file_name() = ".";

  // In general we don't write the partition information.
  
  // If we have p levels and the user has not already provided
  // a file name then write to "."
  if ((n_p_levels > 1) &&
      (this->polynomial_level_file_name() == "n/a"))
    this->polynomial_level_file_name() = ".";

  // write the header
  if (libMesh::processor_id() == 0)
    {
      io.data (this->version());
      
      io.data (n_elem,  "# number of elements");
      io.data (n_nodes, "# number of nodes");
      
      io.data (this->boundary_condition_file_name(), "# boundary condition specification file");
      io.data (this->subdomain_map_file_name(),      "# subdomain id specification file");
      io.data (this->partition_map_file_name(),      "# processor id specification file");
      io.data (this->polynomial_level_file_name(),   "# p-level specification file");      
    }

  // write connectivity
  this->write_serialized_connectivity (io, n_elem);

  // write the nodal locations
  this->write_serialized_nodes (io, n_nodes);

  // write the boundary condition information
  this->write_serialized_bcs (io, n_bcs);

  STOP_LOG("write()","XdrIO");
  
  // pause all processes until the write ends -- this will 
  // protect for the pathological case where a write is 
  // followed immediately by a read.  The write must be
  // guaranteed to complete first.
  io.close();
  Parallel::barrier();
}



void XdrIO::write_serialized_connectivity (Xdr &io, const unsigned int n_elem) const
{
  libmesh_assert (io.writing());
  
  const bool
    write_p_level      = ("." == this->polynomial_level_file_name()),
    write_partitioning = ("." == this->partition_map_file_name()),
    write_subdomain_id = ("." == this->subdomain_map_file_name());

  // convenient reference to our mesh
  const MeshBase &mesh = MeshOutput<MeshBase>::mesh();
  libmesh_assert (n_elem == mesh.n_elem());

  // We will only write active elements and their parents.
  const unsigned int n_active_levels = MeshTools::n_active_levels (mesh);
  std::vector<unsigned int> 
    n_global_elem_at_level(n_active_levels), n_local_elem_at_level(n_active_levels);

  MeshBase::const_element_iterator it  = mesh.local_elements_end(), end=it;

  // Find the number of local and global elements at each level
  for (unsigned int level=0, tot_n_elem=0; level<n_active_levels; level++)
    {
      it  = mesh.local_level_elements_begin(level);
      end = mesh.local_level_elements_end(level);

      n_local_elem_at_level[level] = n_global_elem_at_level[level] = MeshTools::n_elem(it, end);

      Parallel::sum(n_global_elem_at_level[level]);
      tot_n_elem += n_global_elem_at_level[level];
      libmesh_assert (n_global_elem_at_level[level] <= n_elem);
      libmesh_assert (tot_n_elem <= n_elem);
    }

  std::vector<unsigned int> 
    xfer_conn, recv_conn, output_buffer, 
    n_elem_on_proc(libMesh::n_processors()), 
    processor_offsets(libMesh::n_processors()), 
    xfer_buf_sizes(libMesh::n_processors());

  typedef std::map<unsigned int, std::pair<unsigned int, unsigned int> > id_map_type;
  id_map_type parent_id_map, child_id_map;

  unsigned int my_next_elem=0, next_global_elem=0;

  //-------------------------------------------
  // First write the level-0 elements directly.
  it  = mesh.local_level_elements_begin(0);
  end = mesh.local_level_elements_end(0);
  for (; it != end; ++it)
    {
      pack_element (xfer_conn, *it);
      parent_id_map[(*it)->id()] = std::make_pair(libMesh::processor_id(), 
						  my_next_elem++);
    }
  xfer_conn.push_back(my_next_elem); // toss in the number of elements transferred.

  unsigned int my_size = xfer_conn.size();
  Parallel::gather (0, my_next_elem, n_elem_on_proc);
  Parallel::gather (0, my_size,      xfer_buf_sizes);

  processor_offsets[0] = 0;
  for (unsigned int pid=1; pid<libMesh::n_processors(); pid++)
    processor_offsets[pid] = processor_offsets[pid-1] + n_elem_on_proc[pid-1];

  // All processors send their xfer buffers to processor 0.
  Parallel::request request_handle;
  Parallel::isend (0, xfer_conn, request_handle);
  
  // Processor 0 will receive the data and write out the elements. 
  if (libMesh::processor_id() == 0)
    {
      // Write the number of elements at this level.      
      {
	std::string comment = "# n_elem at level 0", legend  = ", [ type ";
	if (write_partitioning)
	  legend += "pid ";
	if (write_subdomain_id)
	  legend += "sid ";
	if (write_p_level)
	  legend += "p_level ";
	legend += "(n0 ... nN-1) ]";
	comment += legend;
	io.data (n_global_elem_at_level[0], comment.c_str());
      }

      for (unsigned int pid=0; pid<libMesh::n_processors(); pid++)
	{
	  recv_conn.resize(xfer_buf_sizes[pid]);
#ifdef HAVE_MPI
	  Parallel::recv (pid, recv_conn);	
#else
	  libmesh_assert (libMesh::n_processors() == 1);
	  libmesh_assert (libMesh::processor_id() == pid);
	  recv_conn = xfer_conn;
#endif
	  // at a minimum, the buffer should contain the number of elements,
	  // which could be 0.
	  libmesh_assert (!recv_conn.empty());
	  
	  const unsigned int n_elem_received = recv_conn.back();
	  std::vector<unsigned int>::const_iterator it = recv_conn.begin();
	  
	  for (unsigned int elem=0; elem<n_elem_received; elem++, next_global_elem++)
	    {
	      output_buffer.clear();
	      const unsigned int n_nodes = *it; ++it;
	      output_buffer.push_back(*it);     /* type       */ ++it;	      
	      /*output_buffer.push_back(*it);*/ /* id         */ ++it;
	      
	      if (write_partitioning)
		output_buffer.push_back(*it); /* processor id */ ++it;

	      if (write_subdomain_id)
		output_buffer.push_back(*it); /* subdomain id */ ++it;

#ifdef ENABLE_AMR
	      if (write_p_level)
		output_buffer.push_back(*it); /* p level      */ ++it;
#endif
	      for (unsigned int node=0; node<n_nodes; node++, ++it)
		output_buffer.push_back(*it);
	      
	      io.data_stream (&output_buffer[0], output_buffer.size(), output_buffer.size());
	    }       	
	}
    }
  Parallel::wait (request_handle);

#ifdef ENABLE_AMR  
  //--------------------------------------------------------------------
  // Next write the remaining elements indirectly through their parents.
  // This will insure that the children are written in the proper order
  // so they can be reconstructed properly.
  for (unsigned int level=1, my_n_elem_written_at_level=0; level<n_active_levels; level++)
    {
      xfer_conn.clear();

      it  = mesh.local_level_elements_begin(level-1);
      end = mesh.local_level_elements_end  (level-1);
      
      for (my_n_elem_written_at_level=0; it != end; ++it)
	if (!(*it)->active()) // we only want the parents elements at this level, and