xdr_io.c

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

	  {                   // there is no direct iterator for this obscure use
	    const Elem *parent = *it;
	    id_map_type::iterator pos = parent_id_map.find(parent->id());
	    libmesh_assert (pos != parent_id_map.end());
	    const unsigned int parent_pid = pos->second.first;
	    const unsigned int parent_id  = pos->second.second;
	    parent_id_map.erase(pos);

	    for (unsigned int c=0; c<parent->n_children(); c++, my_next_elem++)
	      {
		const Elem *child = parent->child(c);
		pack_element (xfer_conn, child, parent_id, parent_pid);

		// this aproach introduces the possibility that we write
		// non-local elements.  These elements may well be parents
		// at the next step
		child_id_map[child->id()] = std::make_pair (child->processor_id(), 
							    my_n_elem_written_at_level++);
	      }
	  }
      xfer_conn.push_back(my_n_elem_written_at_level);
      my_size = xfer_conn.size();
      Parallel::gather (0, my_size,   xfer_buf_sizes);
      Parallel::isend  (0, xfer_conn, request_handle);
      
      // Processor 0 will receive the data and write the elements.
      if (libMesh::processor_id() == 0)
	{
	  // Write the number of elements at this level.      
	  {
	    char buf[80];
	    std::sprintf(buf, "# n_elem at level %d", level);
	    std::string comment(buf), legend  = ", [ type parent ";

	    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[level], 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;
		  
		  const unsigned int parent_local_id = *it; ++it; 
		  const unsigned int parent_pid      = *it; ++it; 
	      	      
		  output_buffer.push_back (parent_local_id+processor_offsets[parent_pid]);

		  if (write_partitioning)
		    output_buffer.push_back(*it); /* processor id */ ++it;

		  if (write_subdomain_id)
		    output_buffer.push_back(*it); /* subdomain id */ ++it;
		  
		  if (write_p_level)
		    output_buffer.push_back(*it); /* p level       */ ++it;
		  
		  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);
  
      // update the processor_offsets
      processor_offsets[0] = processor_offsets.back() + n_elem_on_proc.back();	  
      Parallel::gather (0, my_n_elem_written_at_level, n_elem_on_proc);
      for (unsigned int pid=1; pid<libMesh::n_processors(); pid++)
	processor_offsets[pid] = processor_offsets[pid-1] + n_elem_on_proc[pid-1];	    

      // Now, at the next level we will again iterate over local parents.  However, 
      // those parents may have been written by other processors (at this step), 
      // so we need to gather them into our *_id_maps.
      {
	std::vector<std::vector<unsigned int> > requested_ids(libMesh::n_processors());
	std::vector<unsigned int> request_to_fill;	

	it  = mesh.local_level_elements_begin(level);
	end = mesh.local_level_elements_end(level);

	for (; it!=end; ++it)
	  if (!child_id_map.count((*it)->id()))
	    {
	      libmesh_assert ((*it)->parent()->processor_id() != libMesh::processor_id());
	      requested_ids[(*it)->parent()->processor_id()].push_back((*it)->id());
	    }

	// Next set the child_ids 
	for (unsigned int p=1; p != libMesh::n_processors(); ++p)
	  {
	    // Trade my requests with processor procup and procdown
	    unsigned int procup = (libMesh::processor_id() + p) %
                                   libMesh::n_processors();
	    unsigned int procdown = (libMesh::n_processors() +
                                     libMesh::processor_id() - p) %
	                             libMesh::n_processors();

	    Parallel::send_receive(procup, requested_ids[procup],
				   procdown, request_to_fill);

	    // Fill those requests by overwriting the requested ids
	    for (unsigned int i=0; i<request_to_fill.size(); i++)	   
	      {
		libmesh_assert (child_id_map.count(request_to_fill[i]));
		libmesh_assert (child_id_map[request_to_fill[i]].first == procdown);

		request_to_fill[i] = child_id_map[request_to_fill[i]].second;
	      }

	    // Trade back the results
	    std::vector<unsigned int> filled_request;
	    Parallel::send_receive(procdown, request_to_fill,
				   procup, filled_request);

	    libmesh_assert (filled_request.size() == requested_ids[procup].size());
	    
	    for (unsigned int i=0; i<filled_request.size(); i++)
	      child_id_map[requested_ids[procup][i]] = 
		std::make_pair (procup,
				filled_request[i]);
	  }
	// overwrite the parent_id_map with the child_id_map, but
	// use std::map::swap() for efficiency.
	parent_id_map.swap(child_id_map); /**/ child_id_map.clear();
      }
    }
#endif // ENABLE_AMR
  if (libMesh::processor_id() == 0)
    libmesh_assert (next_global_elem == n_elem);
  
}



void XdrIO::write_serialized_nodes (Xdr &io, const unsigned int n_nodes) const
{
  // convenient reference to our mesh
  const MeshBase &mesh = MeshOutput<MeshBase>::mesh();
  libmesh_assert (n_nodes == mesh.n_nodes());

  std::vector<unsigned int> xfer_ids;
  std::vector<Real>         xfer_coords, &coords=xfer_coords;

  std::vector<std::vector<unsigned int> > recv_ids   (libMesh::n_processors());;
  std::vector<std::vector<Real> >         recv_coords(libMesh::n_processors());

  unsigned int n_written=0;

  for (unsigned int blk=0, last_node=0; last_node<n_nodes; blk++)
    {
      const unsigned int first_node = blk*io_blksize;
                          last_node = std::min((blk+1)*io_blksize, n_nodes);

      // Build up the xfer buffers on each processor
      MeshBase::const_node_iterator
	it  = mesh.local_nodes_begin(),
	end = mesh.local_nodes_end();

      xfer_ids.clear(); xfer_coords.clear();
      
      for (; it!=end; ++it)
	if (((*it)->id() >= first_node) && // node in [first_node, last_node)
	    ((*it)->id() <  last_node))
	  {
	    xfer_ids.push_back((*it)->id());
	    const Point &p = **it;
	    xfer_coords.push_back(p(0));
	    xfer_coords.push_back(p(1));
	    xfer_coords.push_back(p(2));
	  }

      //-------------------------------------
      // Send the xfer buffers to processor 0
      std::vector<unsigned int> ids_size, coords_size;

      const unsigned int my_ids_size = xfer_ids.size();

      // explicitly gather ids_size
      Parallel::gather (0, my_ids_size, ids_size);

      // infer coords_size on processor 0
      if (libMesh::processor_id() == 0)
	{
	  coords_size.reserve(libMesh::n_processors());
	  for (unsigned int p=0; p<ids_size.size(); p++)
	    coords_size.push_back(3*ids_size[p]);
	}			  

      // Note that we will actually send/receive to ourself if we are
      // processor 0, so let's use nonblocking receives.
      std::vector<Parallel::request>
        id_request_handles(libMesh::n_processors()),
        coord_request_handles(libMesh::n_processors());

#ifdef HAVE_MPI
      const unsigned int id_tag=0, coord_tag=1;
      
      // Post the receives -- do this on processor 0 only.
      if (libMesh::processor_id() == 0)
        for (unsigned int pid=0; pid<libMesh::n_processors(); pid++)
          {
            recv_ids[pid].resize(ids_size[pid]);
            recv_coords[pid].resize(coords_size[pid]);
            
            Parallel::irecv (pid, recv_ids[pid],    id_request_handles[pid],    id_tag);
            Parallel::irecv (pid, recv_coords[pid], coord_request_handles[pid], coord_tag);
          }
      
      // Send -- do this on all processors.
      Parallel::send(0, xfer_ids,    id_tag);
      Parallel::send(0, xfer_coords, coord_tag);
#else
      // On one processor there's nothing to send
      recv_ids[0]    = xfer_ids;
      recv_coords[0] = xfer_coords;
#endif

      // -------------------------------------------------------
      // Receive the messages and write the output on processor 0.
      if (libMesh::processor_id() == 0)
        {
          // Wait for all the receives to complete. We have no
          // need for the statuses since we already know the
          // buffer sizes.
	  Parallel::wait (id_request_handles);
	  Parallel::wait (coord_request_handles);
          
          // Write the coordinates in this block.
	  unsigned int tot_id_size=0, tot_coord_size=0;
	  for (unsigned int pid=0; pid<libMesh::n_processors(); pid++)
	    {
	      tot_id_size    += recv_ids[pid].size();
	      tot_coord_size += recv_coords[pid].size();
	    }

	  libmesh_assert (tot_id_size <= std::min(io_blksize, n_nodes));
	  libmesh_assert (tot_coord_size == 3*tot_id_size);

	  coords.resize (tot_coord_size);
	  for (unsigned int pid=0; pid<libMesh::n_processors(); pid++)
	    for (unsigned int idx=0; idx<recv_ids[pid].size(); idx++)
	      {
		const unsigned int local_idx = recv_ids[pid][idx] - first_node;
		libmesh_assert ((3*local_idx+2) < coords.size());
		libmesh_assert ((3*idx+2)      < recv_coords[pid].size());
		
		coords[3*local_idx+0] = recv_coords[pid][3*idx+0];
		coords[3*local_idx+1] = recv_coords[pid][3*idx+1];
		coords[3*local_idx+2] = recv_coords[pid][3*idx+2];

		n_written++;
	      }

	  io.data_stream (coords.empty() ? NULL : &coords[0], coords.size(), 3);
	}
    }
  if (libMesh::processor_id() == 0)
    libmesh_assert (n_written == n_nodes);
}



void XdrIO::write_serialized_bcs (Xdr &io, const unsigned int n_bcs) const
{
  libmesh_assert (io.writing());
  
  if (!n_bcs) return;
  
  // convenient reference to our mesh
  const MeshBase &mesh = MeshOutput<MeshBase>::mesh();
  
  // and our boundary info object
  const BoundaryInfo &boundary_info = *mesh.boundary_info;

  unsigned int n_bcs_out = n_bcs;
  if (libMesh::processor_id() == 0)
    io.data (n_bcs_out, "# number of boundary conditions");
  n_bcs_out = 0;

  std::vector<int> xfer_bcs, recv_bcs;
  std::vector<unsigned int> bc_sizes(libMesh::n_processors());;
  
  // Boundary conditions are only specified for level-0 elements
  MeshBase::const_element_iterator
    it  = mesh.local_level_elements_begin(0),
    end = mesh.local_level_elements_end(0);

  unsigned int n_local_level_0_elem=0;
  for (; it!=end; ++it, n_local_level_0_elem++)
    {
      const Elem *elem = *it;

      for (unsigned int s=0; s<elem->n_sides(); s++)
	if (elem->neighbor(s) == NULL)
	  {
	    const short int bc_id = 
	      boundary_info.boundary_id (elem, s);

	    if (bc_id != BoundaryInfo::invalid_id)
	      {
		xfer_bcs.push_back (n_local_level_0_elem); 
		xfer_bcs.push_back (s) ;
		xfer_bcs.push_back (bc_id);
	      }
	  }
    }

  xfer_bcs.push_back(n_local_level_0_elem);
  unsigned int my_size = xfer_bcs.size();
  Parallel::gather (0, my_size, bc_sizes);
  
  // All processors send their xfer buffers to processor 0
  Parallel::request request_handle;
  Parallel::isend (0, xfer_bcs, request_handle);
  
  // Processor 0 will receive all buffers and write out the bcs
  if (libMesh::processor_id() == 0)
    {
      for (unsigned int pid=0, elem_offset=0; pid<libMesh::n_processors(); pid++)
	{
	  recv_bcs.resize(bc_sizes[pid]);
#ifdef HAVE_MPI
	  Parallel::recv (pid, recv_bcs);
#else
	  libmesh_assert (libMesh::n_processors() == 1);
	  libmesh_assert (libMesh::processor_id() == pid);
	  recv_bcs = xfer_bcs;
#endif
	  const unsigned int n_local_level_0_elem 
	    = recv_bcs.back(); recv_bcs.pop_back();
	  
	  for (unsigned int idx=0; idx<recv_bcs.size(); idx += 3, n_bcs_out++)
	    recv_bcs[idx+0] += elem_offset;

	  io.data_stream (recv_bcs.empty() ? NULL : &recv_bcs[0], recv_bcs.size(), 3);
	  elem_offset += n_local_level_0_elem;
	}    
      libmesh_assert (n_bcs == n_bcs_out);
    }
  Parallel::wait (request_handle);
}