system_io.c

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

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

  // build the ordered nodes and element maps.
  // when writing/reading parallel files we need to iterate
  // over our nodes/elements in order of increasing global id().
  // however, this is not guaranteed to be ordering we obtain
  // by using the node_iterators/element_iterators directly.
  // so build a set, sorted by id(), that provides the ordering.
  // further, for memory economy build the set but then transfer
  // its contents to vectors, which will be sorted.
  std::vector<const DofObject*> ordered_nodes, ordered_elements;
  {    
    std::set<const DofObject*, CompareDofObjectsByID>
      ordered_nodes_set (this->get_mesh().local_nodes_begin(),
			 this->get_mesh().local_nodes_end());
      
      ordered_nodes.insert(ordered_nodes.end(),
			   ordered_nodes_set.begin(),
			   ordered_nodes_set.end());
  }
  {
    std::set<const DofObject*, CompareDofObjectsByID>
      ordered_elements_set (this->get_mesh().local_elements_begin(),
			    this->get_mesh().local_elements_end());
      
      ordered_elements.insert(ordered_elements.end(),
			      ordered_elements_set.begin(),
			      ordered_elements_set.end());
  }
  
  std::vector<Number> io_buffer;
  
  // 9.)
  //
  // Actually read the solution components
  // for the ith system to disk
  io.data(io_buffer);
  
  const unsigned int sys_num = this->number();
  const unsigned int n_vars  = this->n_vars();

  unsigned int cnt=0;
  
  // Loop over each variable and each node, and read out the value.
  for (unsigned int var=0; var<n_vars; var++)
    {
      // First read the node DOF values
      for (std::vector<const DofObject*>::const_iterator
	     it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
	for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
	  {
	    libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
		    DofObject::invalid_id);
	    libmesh_assert (cnt < io_buffer.size());
	    this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
	  }

      // Then read the element DOF values
      for (std::vector<const DofObject*>::const_iterator
	     it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
	for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
	  {
	    libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
		    DofObject::invalid_id);
	    libmesh_assert (cnt < io_buffer.size());
	    this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
	  }      
    }
  
  // Only read additional vectors if wanted  
  if (read_additional_data)
    {	  
      std::map<std::string, NumericVector<Number>* >::const_iterator
	pos = _vectors.begin();
  
      for(; pos != this->_vectors.end(); ++pos)
        {
	  cnt=0;
	  io_buffer.clear();
	  
	  // 10.)
	  //
	  // Actually read the additional vector components
	  // for the ith system to disk
	  io.data(io_buffer);
	  
	  // Loop over each variable and each node, and read out the value.
	  for (unsigned int var=0; var<n_vars; var++)
	    {
	      // First read the node DOF values
	      for (std::vector<const DofObject*>::const_iterator
		     it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
		for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
		  {
		    libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
			    DofObject::invalid_id);
		    libmesh_assert (cnt < io_buffer.size());
		    this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
		  }	     
	      
	      // Then read the element DOF values
	      for (std::vector<const DofObject*>::const_iterator
		     it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
		for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
		  {
		    libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
			    DofObject::invalid_id);
		    libmesh_assert (cnt < io_buffer.size());
		    this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
		  }	      
	    }
	}
    }
}



void System::read_serialized_data (Xdr& io,
				   const bool read_additional_data)
{
  // This method implements the input 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.)
  parallel_only();
  std::string comment;

  // 10.)
  // Read the global solution vector
  {
    this->read_serialized_vector(io, *this->solution); 

    // get the comment
    if (libMesh::processor_id() == 0)
      io.comment (comment);  
  }
  
  // 11.)
  // Only read additional vectors if wanted  
  if (read_additional_data)
    {	  
      std::map<std::string, NumericVector<Number>* >::const_iterator
	pos = _vectors.begin();
  
      for(; pos != this->_vectors.end(); ++pos)
        {
	  this->read_serialized_vector(io, *pos->second);

	  // get the comment
	  if (libMesh::processor_id() == 0)
	    io.comment (comment);	  
	    
	}
    }
}



template <typename iterator_type>
unsigned int System::read_serialized_blocked_dof_objects (const unsigned int var,
							  const unsigned int n_objects,
							  const iterator_type begin,
							  const iterator_type end,
							  Xdr &io,
							  NumericVector<Number> &vec) const
{
  const unsigned int sys_num = this->number();
  
  std::vector<Number> input_buffer;        // buffer to hold the input block read from io.               
  std::vector<Number> local_values;
  std::vector<std::vector<unsigned int> >  // The IDs from each processor which map to the objects
    recv_ids(libMesh::n_processors());     //  read in the current block
  std::vector<unsigned int> idx_map;       // Reordering map to traverse entry-wise rather than processor-wise
  unsigned int n_assigned_vals = 0;        // the number of values assigned, this will be returned.
  
  //-----------------------------------
  // Collect the values for all objects
  unsigned int first_object=0, last_object=0;

  for (unsigned int blk=0; last_object<n_objects; blk++)
    {
      //std::cout << "Reading object block " << blk << std::endl;
      
      // Each processor should build up its transfer buffers for its
      // local objects in [first_object,last_object).
      first_object = blk*io_blksize;
      last_object  = std::min((blk+1)*io_blksize,n_objects);
      
      // Clear the transfer buffers for this block.
      recv_ids[libMesh::processor_id()].clear();
      unsigned int n_local_dofs=0;
      for (iterator_type it=begin; it!=end; ++it)
	if (((*it)->id() >= first_object) && // object in [first_object,last_object)
	    ((*it)->id() <   last_object) &&
	    (*it)->n_comp(sys_num,var))      // var has a nonzero # of components on this object
	  {
	    recv_ids[libMesh::processor_id()].push_back((*it)->id());
	    recv_ids[libMesh::processor_id()].push_back((*it)->n_comp(sys_num, var));
	    recv_ids[libMesh::processor_id()].push_back(n_local_dofs);
	    n_local_dofs += (*it)->n_comp(sys_num, var);	    
	  }
      
      // Get the recv_ids for all other processors.
      {
	const unsigned int curr_vec_size = recv_ids[libMesh::processor_id()].size();
	std::vector<unsigned int> recv_id_sizes(libMesh::n_processors());
	Parallel::allgather(curr_vec_size, recv_id_sizes);
	for (unsigned int pid=0; pid<libMesh::n_processors(); pid++)
	  {
	    recv_ids[pid].resize(recv_id_sizes[pid]);
	    Parallel::broadcast(recv_ids[pid], pid);
	  }
      }

      // create the idx map for all processors -- this will match the ordering
      // in the input buffer chunk which we are about to read.
      idx_map.resize(3*io_blksize); std::fill (idx_map.begin(), idx_map.end(), libMesh::invalid_uint);
      unsigned int tot_n_comp=0;
      for (unsigned int pid=0; pid<libMesh::n_processors(); pid++)
	for (unsigned int idx=0; idx<recv_ids[pid].size(); idx+=3)
	  {
	    const unsigned int local_idx = recv_ids[pid][idx+0]-first_object;
	    libmesh_assert (local_idx < std::min(io_blksize,n_objects));
	    const unsigned int n_comp    = recv_ids[pid][idx+1];
	    const unsigned int start_pos = recv_ids[pid][idx+2];
	    
	    idx_map[3*local_idx+0] = pid;
	    idx_map[3*local_idx+1] = n_comp;
	    idx_map[3*local_idx+2] = start_pos; // this tells us where the values
	                                        // for this object will live in the local_values buffer
	    tot_n_comp += n_comp;
	  }
      
      // Processor 0 will read the block from the buffer stream and send it to all other processors
      input_buffer.resize (tot_n_comp);
      if (libMesh::processor_id() == 0) io.data_stream(input_buffer.empty() ? NULL : &input_buffer[0], tot_n_comp);
      Parallel::broadcast(input_buffer);

      // Extract the values corresponding to the local objects from the input_buffer
      // and place them into the local_values temporary buffer.
      {
	unsigned int processed_size=0;
	std::vector<Number>::const_iterator next_value = input_buffer.begin();
	local_values.clear(); local_values.resize(n_local_dofs);

	for (unsigned int idx=0; idx<idx_map.size(); idx+=3)
	  if (idx_map[idx] != libMesh::invalid_uint) // this could happen when an object
	    {                                        // has no components for the current variable
	      const unsigned int pid       = idx_map[idx+0];
	      const unsigned int n_comp    = idx_map[idx+1];
	      const unsigned int start_pos = idx_map[idx+2];
	      
	      for (unsigned int comp=0; comp<n_comp; comp++)
		{
		  libmesh_assert (next_value != input_buffer.end());
		  if (pid == libMesh::processor_id())
		    {
		      libmesh_assert ((start_pos+comp) < local_values.size());
		      local_values[start_pos+comp] = *next_value;
		    }
		  ++next_value;
		  ++processed_size;				
		}
	    }
	libmesh_assert (processed_size == input_buffer.size());
      }
      
      // A subset of the components (potentially null set) will match our objects in
      // [first_object,last_object), and we will assign the corresponding values from
      // the local_values buffer.
      for (iterator_type it=begin; it!=end; ++it)
	if (((*it)->id() >= first_object) && // object in [first_object,last_object)
	    ((*it)->id() <   last_object) &&
	    (*it)->n_comp(sys_num,var))      // var has a nonzero # of components on this object
	  {
	    const unsigned int local_idx = (*it)->id()-first_object;
	    libmesh_assert (local_idx < std::min(io_blksize,n_objects));

#ifndef NDEBUG
	    // We only need to check the pid when asserts are active
	    const unsigned int pid       = idx_map[3*local_idx+0];
#endif
	    libmesh_assert (pid == libMesh::processor_id());
	    
	    const unsigned int n_comp    = idx_map[3*local_idx+1];
	    const unsigned int start_pos = idx_map[3*local_idx+2];
	    
	    libmesh_assert (n_comp == (*it)->n_comp(sys_num, var));
	    
	    for (unsigned int comp=0; comp<n_comp; comp++)
	      {
		libmesh_assert ((start_pos+comp) < local_values.size());
		const Number &value = local_values[start_pos+comp];				
		const unsigned int dof_index = (*it)->dof_number (sys_num, var, comp);
		libmesh_assert (dof_index >= vec.first_local_index());
		libmesh_assert (dof_index <  vec.last_local_index());
		//std::cout << "di=" << dof_index << ", val=" << value << std::endl;
		vec.set (dof_index, value);
		++n_assigned_vals;
	      }
	  }
    }
      
  return n_assigned_vals;
}



void System::read_serialized_vector (Xdr& io, NumericVector<Number>& vec)
{
  parallel_only();

#ifndef NDEBUG
  // In parallel we better be reading a parallel vector -- if not
  // we will not set all of its components below!!
  if (libMesh::n_processors() > 1)
    libmesh_assert (vec.size() != vec.local_size());
   
  // If this is not the same on all processors we're in trouble!
  Parallel::verify(io_blksize);
#endif
  
  libmesh_assert (io.reading());

  // vector length
  unsigned int vector_length=0, n_assigned_vals=0;

  // Get the buffer size
  if (libMesh::processor_id() == 0) io.data(vector_length, "# vector length");
  Parallel::broadcast(vector_length);
  
  // Loop over each variable in the system, and then each node/element in the mesh.
  for (unsigned int var=0; var<this->n_vars(); var++)
    {      
      //---------------------------------
      // Collect the values for all nodes
      n_assigned_vals +=
	this->read_serialized_blocked_dof_objects (var,
						   this->get_mesh().n_nodes(),
						   this->get_mesh().local_nodes_begin(),
						   this->get_mesh().local_nodes_end(),
						   io,
						   vec);
      
      
      //------------------------------------
      // Collect the values for all elements
      n_assigned_vals +=
	this->read_serialized_blocked_dof_objects (var,
						   this->get_mesh().n_elem(),
						   this->get_mesh().local_elements_begin(),
						   this->get_mesh().local_elements_end(),
						   io,
						   vec);
    } // end variable loop

  vec.close();

#ifdef DEBUG
  Parallel::sum (n_assigned_vals);
#endif
  libmesh_assert (n_assigned_vals == vector_length);
}



void System::write_header (Xdr& io,
			   const bool write_additional_data) const
{
  /**
   * This method implements the output of a
   * System object, embedded in the output of
   * an EquationSystems<T_sys>.  This warrants some