xdr_io.c

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

void XdrIO::read (const std::string& name)
{
  // Only open the file on processor 0 -- this is especially important because
  // there may be an underlying bzip/bunzip going on, and multiple simultaneous
  // calls will produce a race condition.
  Xdr io (libMesh::processor_id() == 0 ? name : "", this->binary() ? DECODE : READ);
   
  // convenient reference to our mesh
  MeshBase &mesh = MeshInput<MeshBase>::mesh();

  // get the version string.   
  if (libMesh::processor_id() == 0)    
    io.data (this->version());
  Parallel::broadcast (this->version());  
  
  // note that for "legacy" files the first entry is an
  // integer -- not a string at all.
  this->legacy() = !(this->version().find("libMesh") < this->version().size());

  // Check for a legacy version format.
  if (this->legacy())
    {
      io.close();
      LegacyXdrIO(mesh, this->binary()).read(name);
      return;
    }

  START_LOG("read()","XdrIO");
  
  unsigned int n_elem, n_nodes;
  if (libMesh::processor_id() == 0)
    {
      io.data (n_elem);
      io.data (n_nodes);
      io.data (this->boundary_condition_file_name()); // std::cout << "bc_file="  << this->boundary_condition_file_name() << std::endl;
      io.data (this->subdomain_map_file_name());      // std::cout << "sid_file=" << this->subdomain_map_file_name()      << std::endl;
      io.data (this->partition_map_file_name());      // std::cout << "pid_file=" << this->partition_map_file_name()      << std::endl;
      io.data (this->polynomial_level_file_name());   // std::cout << "pl_file="  << this->polynomial_level_file_name()   << std::endl;
    }

  //TODO:[BSK] a little extra effort here could change this to two broadcasts...
  Parallel::broadcast (n_elem);
  Parallel::broadcast (n_nodes);
  Parallel::broadcast (this->boundary_condition_file_name());
  Parallel::broadcast (this->subdomain_map_file_name());
  Parallel::broadcast (this->partition_map_file_name());
  Parallel::broadcast (this->polynomial_level_file_name());
  
  // Tell the mesh how many nodes/elements to expect. Depending on the mesh type,
  // this may allow for efficient adding of nodes/elements.
  mesh.reserve_elem(n_elem);
  mesh.reserve_nodes(n_nodes);

  // read connectivity
  this->read_serialized_connectivity (io, n_elem);

  // read the nodal locations
  this->read_serialized_nodes (io, n_nodes);  

  // read the boundary conditions
  this->read_serialized_bcs (io);

  STOP_LOG("read()","XdrIO");
   
  // set the node processor ids
  Partitioner::set_node_processor_ids(mesh);
}



void XdrIO::read_serialized_connectivity (Xdr &io, const unsigned int n_elem)
{
  libmesh_assert (io.reading());

  if (!n_elem) return;
  
  const bool
    read_p_level      = ("." == this->polynomial_level_file_name()),
    read_partitioning = ("." == this->partition_map_file_name()),
    read_subdomain_id = ("." == this->subdomain_map_file_name());
  
  // convenient reference to our mesh
  MeshBase &mesh = MeshInput<MeshBase>::mesh();

  std::vector<unsigned int> conn, input_buffer(100 /* oversized ! */);

  int level=-1;
  
  for (unsigned int blk=0, first_elem=0, last_elem=0, n_elem_at_level=0, n_processed_at_level=0;
       last_elem<n_elem; blk++)
    {
      first_elem = blk*io_blksize;
      last_elem  = std::min((blk+1)*io_blksize, n_elem);

      conn.clear();

      if (libMesh::processor_id() == 0)
	for (unsigned int e=first_elem; e<last_elem; e++, n_processed_at_level++)
	  {
	    if (n_processed_at_level == n_elem_at_level)
	      {
		// get the number of elements to read at this level
		io.data (n_elem_at_level);
		n_processed_at_level = 0;
		level++;
	      }
	    
	    // get the element type, 
	    io.data_stream (&input_buffer[0], 1);

	    // maybe the parent
	    if (level)
	      io.data_stream (&input_buffer[1], 1);
	    else
	      input_buffer[1] = libMesh::invalid_uint;

	    // maybe the processor id
	    if (read_partitioning)
	      io.data_stream (&input_buffer[2], 1);
	    else
	      input_buffer[2] = 0;

	    // maybe the subdomain id
	    if (read_subdomain_id)
	      io.data_stream (&input_buffer[3], 1);
	    else 
	      input_buffer[3] = 0;

	    // maybe the p level
	    if (read_p_level)
	      io.data_stream (&input_buffer[4], 1);
	    else
	      input_buffer[4] = 0;
	    
	    // and all the nodes
	    libmesh_assert (5+Elem::type_to_n_nodes_map[input_buffer[0]] < input_buffer.size());
	    io.data_stream (&input_buffer[5], Elem::type_to_n_nodes_map[input_buffer[0]]);
	    conn.insert (conn.end(),
			 input_buffer.begin(),
			 input_buffer.begin() + 5 + Elem::type_to_n_nodes_map[input_buffer[0]]);
	  }
      
      unsigned int conn_size = conn.size();
      Parallel::broadcast(conn_size);	    
      conn.resize (conn_size);
      Parallel::broadcast (conn);

      // All processors now have the connectivity for this block.
      std::vector<unsigned int>::const_iterator it = conn.begin();      
      for (unsigned int e=first_elem; e<last_elem; e++)
	{
	  const ElemType elem_type        = static_cast<ElemType>(*it); ++it;
	  const unsigned int parent_id    = *it; ++it;
	  const unsigned int processor_id = *it; ++it;
	  const unsigned int subdomain_id = *it; ++it;
#ifdef ENABLE_AMR
	  const unsigned int p_level      = *it;
#endif
	  ++it;

	  Elem *parent = (parent_id == libMesh::invalid_uint) ? NULL : mesh.elem(parent_id);

	  Elem *elem = Elem::build (elem_type, parent).release();
	  elem->set_id() = e;
	  elem->processor_id() = processor_id;
	  elem->subdomain_id() = subdomain_id;
#ifdef ENABLE_AMR
	  elem->hack_p_level(p_level);

	  if (parent)
	    {
	      parent->add_child(elem);
	      parent->set_refinement_flag (Elem::INACTIVE);
	      elem->set_refinement_flag   (Elem::JUST_REFINED);
	    }
#endif

	  for (unsigned int n=0; n<elem->n_nodes(); n++, ++it)
	    {
	      const unsigned int global_node_number = *it;

	      elem->set_node(n) = 
		mesh.add_point (Point(), global_node_number);
	    }
	  
	  mesh.add_elem(elem);
	}
    }
}



void XdrIO::read_serialized_nodes (Xdr &io, const unsigned int n_nodes)
{
  libmesh_assert (io.reading());
  
  // convenient reference to our mesh
  MeshBase &mesh = MeshInput<MeshBase>::mesh();
  
  if (!mesh.n_nodes()) return;

  // At this point the elements have been read from file and placeholder nodes
  // have been assigned.  These nodes, however, do not have the proper (x,y,z)
  // locations.  This method will read all the nodes from disk, and each processor
  // can then grab the individual values it needs.

  // build up a list of the nodes contained in our local mesh.  These are the nodes
  // stored on the local processor whose (x,y,z) values need to be corrected.
  std::vector<unsigned int> needed_nodes; needed_nodes.reserve (mesh.n_nodes());
  {
    MeshBase::node_iterator
      it  = mesh.nodes_begin(),
      end = mesh.nodes_end();

    for (; it!=end; ++it)
      needed_nodes.push_back((*it)->id());

    std::sort (needed_nodes.begin(), needed_nodes.end());

    // We should not have any duplicate node->id()s
    libmesh_assert (std::unique(needed_nodes.begin(), needed_nodes.end()) == needed_nodes.end());
  }

  // Get the nodes in blocks.
  std::vector<Real> coords;
  std::pair<std::vector<unsigned int>::iterator,
            std::vector<unsigned int>::iterator> pos;
  pos.first = needed_nodes.begin();

  for (unsigned int blk=0, first_node=0, last_node=0; last_node<n_nodes; blk++)
    {
      first_node = blk*io_blksize;
      last_node  = std::min((blk+1)*io_blksize, n_nodes);
      
      coords.resize(3*(last_node - first_node));
	
      if (libMesh::processor_id() == 0)
	io.data_stream (coords.empty() ? NULL : &coords[0], coords.size());

      // For large numbers of processors the majority of processors at any given
      // block may not actually need these data.  It may be worth profiling this,
      // although it is expected that disk IO will be the bottleneck
      Parallel::broadcast (coords);
      
      for (unsigned int n=first_node, idx=0; n<last_node; n++, idx+=3)
	{
	  // first see if we need this node.  use pos.first as a smart lower
	  // bound, this will ensure that the size of the searched range
	  // decreases as we match nodes.
	  pos = std::equal_range (pos.first, needed_nodes.end(), n);

	  if (pos.first != pos.second) // we need this node.
	    {
	      libmesh_assert (*pos.first == n);
	      mesh.node(n) = 
		Point (coords[idx+0],
		       coords[idx+1],
		       coords[idx+2]);

	    }
	}	  
    }    
}



void XdrIO::read_serialized_bcs (Xdr &io)
{
  if (this->boundary_condition_file_name() == "n/a") return;

  libmesh_assert (io.reading());
  
  // convenient reference to our mesh
  MeshBase &mesh = MeshInput<MeshBase>::mesh();
  
  // and our boundary info object
  BoundaryInfo &boundary_info = *mesh.boundary_info;

  std::vector<ElemBCData> elem_bc_data;
  std::vector<int> input_buffer;

  unsigned int n_bcs=0;
  if (libMesh::processor_id() == 0) 
    io.data (n_bcs);
  Parallel::broadcast (n_bcs);

  for (unsigned int blk=0, first_bc=0, last_bc=0; last_bc<n_bcs; blk++)
    {
      first_bc = blk*io_blksize;
      last_bc  = std::min((blk+1)*io_blksize, n_bcs);

      input_buffer.resize (3*(last_bc - first_bc));

      if (libMesh::processor_id() == 0)
	io.data_stream (input_buffer.empty() ? NULL : &input_buffer[0], input_buffer.size());
      
      Parallel::broadcast (input_buffer);
      elem_bc_data.clear(); /**/ elem_bc_data.reserve (input_buffer.size()/3);

      // convert the input_buffer to ElemBCData to facilitate searching
      for (unsigned int idx=0; idx<input_buffer.size(); idx+=3)
	elem_bc_data.push_back (ElemBCData(input_buffer[idx+0],
					   input_buffer[idx+1],
					   input_buffer[idx+2]));
      input_buffer.clear();
      // note that while the files *we* write should already be sorted by
      // element id this is not necessarily guaranteed.
      std::sort (elem_bc_data.begin(), elem_bc_data.end());

      MeshBase::const_element_iterator
	it  = mesh.level_elements_begin(0),
	end = mesh.level_elements_end(0);

      // Look for BCs in this block for all the level-0 elements we have 
      // (not just local ones).  Do this by finding all the entries
      // in elem_bc_data whose elem_id match the ID of the current element.
      // We cannot rely on NULL neighbors at this point since the neighbor
      // data structure has not been initialized.
      for (std::pair<std::vector<ElemBCData>::iterator,
	             std::vector<ElemBCData>::iterator> pos; it!=end; ++it)
#if defined(__SUNPRO_CC) || defined(__PGI)
	for (pos = std::equal_range (elem_bc_data.begin(), elem_bc_data.end(), (*it)->id(), CompareIntElemBCData());
#else
	for (pos = std::equal_range (elem_bc_data.begin(), elem_bc_data.end(), (*it)->id());
#endif
	     pos.first != pos.second; ++pos.first)
	  {
	    libmesh_assert (pos.first->elem_id == (*it)->id());
	    libmesh_assert (pos.first->side < (*it)->n_sides());

	    boundary_info.add_side (*it, pos.first->side, pos.first->bc_id);
	  }	
    }
}



void XdrIO::pack_element (std::vector<unsigned int> &conn, const Elem *elem, 
			  const unsigned int parent_id, const unsigned int parent_pid) const
{
  libmesh_assert (elem != NULL);
  libmesh_assert (elem->n_nodes() == Elem::type_to_n_nodes_map[elem->type()]);
  
  conn.push_back(elem->n_nodes());

  conn.push_back (elem->type());
  conn.push_back (elem->id());
  
  if (parent_id != libMesh::invalid_uint)
    {
      conn.push_back (parent_id);
      libmesh_assert (parent_pid != libMesh::invalid_uint);
      conn.push_back (parent_pid);
    }
  
  conn.push_back (elem->processor_id());
  conn.push_back (elem->subdomain_id());
  
#ifdef ENABLE_AMR
  conn.push_back (elem->p_level());
#endif
  
  for (unsigned int n=0; n<elem->n_nodes(); n++)
    conn.push_back (elem->node(n));                 
}