mesh_communication.c

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

	  recv_n_nodes_and_elem_per_proc[5*source_pid+0];
	libmesh_assert (n_nodes_received <= received_nodes.size());
	libmesh_assert (status.size() == n_nodes_received);
	libmesh_assert (recv_node_pair[source_pid]);
	
	for (unsigned int n=0; n<n_nodes_received; n++)
	  {
	    Node *node = received_nodes[n].build_node().release();
	    mesh.insert_node(node); // insert_node works even if the
	  }                         // node already exists in the mesh,
      }                             // in which case it overwrites (x,y,z)
  }

  // Receive BCs for those nodes.
  {
    std::vector<int> received_node_bcs(max_node_bcs_received);

    // again, receive all the messages in whatever order they come.
    for (unsigned int node_comm_step=0; node_comm_step<n_recv_node_bc_pairs; node_comm_step++)
      {
	Parallel::Status status =
	  Parallel::recv (Parallel::any_source,
			  received_node_bcs,
			  /* tag = */ 2);

	libmesh_assert (status.size() == recv_n_nodes_and_elem_per_proc[5*status.source()+3]);

	const unsigned int buffer_size = status.size();
	unsigned int cnt=0;

	while (cnt < buffer_size)
	  {
	    const unsigned int node_id = received_node_bcs[cnt++];
	    const int bc_id            = received_node_bcs[cnt++];

	    libmesh_assert (mesh.node_ptr(node_id));
	    mesh.boundary_info->add_node (mesh.node_ptr(node_id), bc_id);
	  }	
      }
  }


  
  // Receive elements.  Size this array for the largest message.
  {
    std::vector<int> received_elements(max_conn_size_received);

    // Similarly we know how many processors are sending us elements, 
    // but we don't really care in what order we receive them.
    for (unsigned int elem_comm_step=0; elem_comm_step<n_recv_elem_pairs; elem_comm_step++)
      {
	Parallel::Status status =
	  Parallel::recv (Parallel::any_source,
			  received_elements,
			  /* tag = */ 1);

#ifndef NDEBUG
	// Avoid declaring these variables unless asserts are enabled.
	const unsigned int source_pid = status.source();
	const unsigned int n_elem_received = 
	  recv_n_nodes_and_elem_per_proc[5*source_pid+1];
#endif
	libmesh_assert (recv_elem_pair[source_pid]);
	libmesh_assert (recv_n_nodes_and_elem_per_proc[5*source_pid+2] 
		<= received_elements.size());
	libmesh_assert (recv_n_nodes_and_elem_per_proc[5*source_pid+2] 
		== status.size());

	// iterate through the input buffer and add the elements
	const unsigned int xfer_buffer_size = status.size();
	unsigned int cnt=0;
	unsigned int current_elem=0;
	while (cnt < xfer_buffer_size)
	  {
	    // Unpack the element 
	    Elem::PackedElem packed_elem (received_elements.begin()+cnt);
	    
	    // The ParallelMesh::elem(i) member will return NULL if the element
	    // is not in the mesh.  We rely on that here, so it better not change!
	    Elem *elem = mesh.elem(packed_elem.id());
	    
	    // if we already have this element, make sure its properties match
	    // but then go on
	    if (elem)
	      {
		libmesh_assert (elem->level()             == packed_elem.level());
		libmesh_assert (elem->id()                == packed_elem.id());
		libmesh_assert (elem->processor_id()      == packed_elem.processor_id());
		libmesh_assert (elem->subdomain_id()      == packed_elem.subdomain_id());
		libmesh_assert (elem->type()              == packed_elem.type());
#ifdef ENABLE_AMR
		libmesh_assert (elem->p_level()           == packed_elem.p_level());
		libmesh_assert (elem->refinement_flag()   == packed_elem.refinement_flag());
		libmesh_assert (elem->p_refinement_flag() == packed_elem.p_refinement_flag());
		
		if (elem->level() > 0)
		  {
		    libmesh_assert (elem->parent()->id() == static_cast<unsigned int>(packed_elem.parent_id()));
		    libmesh_assert (elem->parent()->child(packed_elem.which_child_am_i()) == elem);
		  }	      
#endif
		libmesh_assert (elem->n_nodes() == packed_elem.n_nodes());
	      }
	    else
	      {
		// We need to add the element.
#ifdef ENABLE_AMR
		// maybe find the parent
		if (packed_elem.level() > 0)
		  {
		    Elem *parent = mesh.elem(packed_elem.parent_id());
		    
		    // Note that we were very careful to construct the send connectivity
		    // so that parents are encountered before children.  So if we get here
		    // and can't find the parent that is a fatal error.
		    if (parent == NULL)
		      {
			std::cerr << "Parent element with ID " << packed_elem.parent_id()
				  << " not found." << std::endl; 
			libmesh_error();
		      }
		    
		    elem = packed_elem.unpack (mesh, parent);
		  }
		else
		  {
		    libmesh_assert (packed_elem.parent_id() == -1);
#endif // ENABLE_AMR
		    elem = packed_elem.unpack (mesh);
#ifdef ENABLE_AMR
		  }
#endif		
		// Good to go.  Add to the mesh.
		libmesh_assert (elem);
		libmesh_assert (elem->n_nodes() == packed_elem.n_nodes());
		mesh.insert_elem(elem);
	      }
	    
	    // properly position cnt for the next element 
	    cnt += Elem::PackedElem::header_size + packed_elem.n_nodes();
	    current_elem++;
	  }
	libmesh_assert (current_elem == n_elem_received);
      }
  }

  // Receive boundary conditions for those elements
  {
    std::vector<int> received_element_bcs(max_elem_bcs_received);
    
    // again, receive all the messages in whatever order they come.
    for (unsigned int elem_comm_step=0; elem_comm_step<n_recv_elem_bc_pairs; elem_comm_step++)
      {
	Parallel::Status status =
	  Parallel::recv (Parallel::any_source,
			  received_element_bcs,
			  /* tag = */ 3);

	libmesh_assert (status.size() == recv_n_nodes_and_elem_per_proc[5*status.source()+4]);

	const unsigned int buffer_size = status.size();
	unsigned int cnt=0;

	while (cnt < buffer_size)
	  {
	    const unsigned int elem_id = received_element_bcs[cnt++];
	    const unsigned int side    = received_element_bcs[cnt++];
	    const int bc_id            = received_element_bcs[cnt++];

	    libmesh_assert (mesh.elem(elem_id));
	    mesh.boundary_info->add_side (mesh.elem(elem_id), side, bc_id);
	  }
      }    
  }

  // Wait for all sends to complete
  Parallel::wait (node_send_requests);
  Parallel::wait (node_bc_requests);
  Parallel::wait (element_send_requests);
  Parallel::wait (element_bc_requests);
  
  // unregister MPI datatypes
  MPI_Type_free (&packed_node_datatype);
  
  STOP_LOG("redistribute()","MeshCommunication");  
}
#endif // HAVE_MPI




void MeshCommunication::broadcast (MeshBase& mesh) const
{
  // Don't need to do anything if there is
  // only one processor.
  if (libMesh::n_processors() == 1)
    return;
  
  this->broadcast_mesh (mesh);
  this->broadcast_bcs  (mesh, *(mesh.boundary_info));
}



#ifdef HAVE_MPI
void MeshCommunication::broadcast_mesh (MeshBase& mesh) const
#else // avoid spurious gcc warnings
void MeshCommunication::broadcast_mesh (MeshBase&) const
#endif
{
  // Don't need to do anything if there is
  // only one processor.
  if (libMesh::n_processors() == 1)
    return;
  
#ifdef HAVE_MPI

  // This function must be run on all processors at once
  parallel_only();

  START_LOG("broadcast_mesh()","MeshCommunication");

  // Explicitly clear the mesh on all but processor 0.
  if (libMesh::processor_id() != 0)
    mesh.clear();
  
  // Get important sizes
  unsigned int n_nodes      = mesh.n_nodes();
  unsigned int n_elem       = mesh.n_elem();
  unsigned int n_levels     = MeshTools::n_levels(mesh);
  unsigned int total_weight = MeshTools::total_weight(mesh);

  // Broadcast the sizes
  {
    std::vector<unsigned int> buf (3);
    
    buf[0] = n_nodes;
    buf[1] = n_elem;
    buf[2] = total_weight;
    
    // Broadcast
    Parallel::broadcast (buf);

    if (libMesh::processor_id() != 0)
      {
	n_nodes      = buf[0];
	n_elem       = buf[1];
	total_weight = buf[2];
      }	
  }  

  // First build up the pts vector which contains
  // the spatial locations of all the nodes      
  {
    std::vector<Real> pts;
	
    // If we are processor 0, we must populate this vector and
    // broadcast it to the other processors.
    if (libMesh::processor_id() == 0)
      {
	pts.reserve (3*n_nodes);

	MeshBase::node_iterator       it     = mesh.nodes_begin();
	const MeshBase::node_iterator it_end = mesh.nodes_end();

	for (; it != it_end; ++it)
	  {
	    libmesh_assert (*it != NULL);
            libmesh_assert (!(*it)->valid_processor_id());
            libmesh_assert ((*it)->id()*3 == pts.size());
	    
	    const Point& p = **it;
	    
	    pts.push_back ( p(0) ); // x
	    pts.push_back ( p(1) ); // y
	    pts.push_back ( p(2) ); // z	  
	  }
      }
    else
      pts.resize (3*n_nodes);

    // Sanity check for all processors
    libmesh_assert (pts.size() == (3*n_nodes));
    
    // Broadcast the pts vector
    Parallel::broadcast (pts);

    // Add the nodes we just received if we are not
    // processor 0.
    if (libMesh::processor_id() != 0)
      {
	libmesh_assert (mesh.n_nodes() == 0);
	
	for (unsigned int i=0; i<pts.size(); i += 3)
	  mesh.add_point (Point(pts[i+0],
				pts[i+1],
				pts[i+2]),
			  i/3);
      }
    
    libmesh_assert (mesh.n_nodes() == n_nodes);
  } // Done distributing the nodes

  
  // Now build up the elements vector which
  // contains the element types and connectivity
  {
    // The conn array contains the information needed to construct each element.
    // Pack all this information into one communication to avoid two latency hits
    // For each element it is of the form
    // [ level p_level r_flag p_flag etype subdomain_id 
    //   self_ID parent_ID which_child node_0 node_1 ... node_n]
    // We cannot use unsigned int because parent_ID can be negative
    std::vector<int> conn;

    // If we are processor 0, we must populate this vector and
    // broadcast it to the other processors.
    if (libMesh::processor_id() == 0)
      {
	conn.reserve (Elem::PackedElem::header_size*n_elem + total_weight);
	
	// We start from level 0. This is a bit simpler than in xdr_io.C
	// because we do not have to worry about economizing by group elements
	// of the same type. Element type is simply specified as an
	// entry in the connectivity vector, "conn".
	// By filling conn in order of levels, parents should exist before children
	// are built when we reconstruct the elements on the other processors.
	
	for (unsigned int level=0; level<=n_levels; ++level)
	  {
	    MeshBase::element_iterator it = mesh.level_elements_begin(level);
	    const MeshBase::element_iterator it_end = mesh.level_elements_end(level);
	    
	    for (; it != it_end; ++it)
	      {
                libmesh_assert (*it);
                libmesh_assert (!(*it)->valid_processor_id());
		const Elem* elem = *it;
		Elem::PackedElem::pack (conn, elem);
	      }
	  }
      }
    else
      conn.resize (Elem::PackedElem::header_size*n_elem + total_weight);
    
    // Sanity check for all processors
    libmesh_assert (conn.size() == (Elem::PackedElem::header_size*n_elem + total_weight));
    
    // Broadcast the element connectivity
    Parallel::broadcast (conn);

    // Build the elements we just received if we are not
    // processor 0.
    if (libMesh::processor_id() != 0)
      {
	libmesh_assert (mesh.n_elem() == 0);
	
	unsigned int cnt = 0;

        // This map keeps track of elements we've previously added to the mesh 
        // to avoid O(n) lookup times for parent pointers.
        std::map<unsigned int, Elem*> parents;

	while (cnt < conn.size())
	  {
	    // Unpack the element
	    Elem::PackedElem packed_elem (conn.begin()+cnt);

	    // Declare the element that we will add
            Elem* elem = NULL;

// 	    // Unpack the element header
// #ifdef ENABLE_AMR
// 	    const int level             = conn[cnt++];
// 	    const int p_level           = conn[cnt++];
// 	    const Elem::RefinementState refinement_flag =
//               static_cast<Elem::RefinementState>(conn[cnt++]);
// 	    const Elem::RefinementState p_refinement_flag =
//               static_cast<Elem::RefinementState>(conn[cnt++]);
// #endif
//             const ElemType elem_type    = static_cast<ElemType>(conn[cnt++]);
// 	    const unsigned int elem_PID = conn[cnt++];
// 	    const int subdomain_ID      = conn[cnt++];
//             const int self_ID           = conn[cnt++];
// #ifdef ENABLE_AMR
//             const int parent_ID         = conn[cnt++];
//             const int which_child       = conn[cnt++];

#ifdef ENABLE_AMR

            if (packed_elem.parent_id() != -1) // Do a log(n) search for the parent
	      {
		Elem* my_parent = 
		  parents.count(packed_elem.parent_id()) ? parents[packed_elem.parent_id()] : NULL;
                
                // If the parent was not previously added, we cannot continue.
                if (my_parent == NULL)
		  {
		    std::cerr << "Parent element with ID " << packed_elem.parent_id()
			      << " not found." << std::endl; 
		    libmesh_error();
		  }
		
                libmesh_assert (my_parent->refinement_flag() == Elem::INACTIVE);		
		
		elem = packed_elem.unpack (mesh, my_parent);
	      }
	    
            else // level 0 element has no parent
	      {
		libmesh_assert (packed_elem.level() == 0);
#endif 		
		elem = packed_elem.unpack (mesh);
#ifdef ENABLE_AMR
	      }
#endif	    
            // Add elem to the map of parents, since it may have
            // children to be added later
            parents.insert(std::make_pair(elem->id(),elem));

	    cnt += Elem::PackedElem::header_size + packed_elem.n_nodes();
	  } // end while cnt < conn.size

        // Iterate in ascending elem ID order
        for (std::map<unsigned int, Elem *>::iterator i =
             parents.begin();
             i != parents.end(); ++i)
          {
            Elem *elem = i->second;
            if (elem)
              mesh.add_elem(elem);
            else
              // We can probably handle this, but we don't expect it
              libmesh_error();
          }

      } // end if iam != cpu 0
    
    
    libmesh_assert (mesh.n_elem() == n_elem);
  } // Done distributing the elements


  STOP_LOG("broadcast_mesh()","MeshCommunication");
  
#else

  // no MPI but multiple processors? Huh??
  libmesh_error();
  
#endif
}



#ifdef HAVE_MPI
void MeshCommunication::broadcast_bcs (const MeshBase& mesh,
				       BoundaryInfo& boundary_info) const
#else // avoid spurious gcc warnings
void MeshCommunication::broadcast_bcs (const MeshBase&,
				       BoundaryInfo&) const
#endif
{
  // Don't need to do anything if there is
  // only one processor.
  if (libMesh::n_processors() == 1)
    return;