parallel_mesh.c

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

  // Clear our elements and nodes
  {
    elem_iterator_imp        it = _elements.begin();
    const elem_iterator_imp end = _elements.end();

    // There is no need to remove the elements from
    // the BoundaryInfo data structure since we
    // already cleared it.
    for (; it != end; ++it)
      delete *it;

    _elements.clear();
  }

  // clear the nodes data structure
  {
    node_iterator_imp it  = _nodes.begin();
    node_iterator_imp end = _nodes.end();

    // There is no need to remove the nodes from
    // the BoundaryInfo data structure since we
    // already cleared it.
    for (; it != end; ++it)
      delete *it;
    
    _nodes.clear();
  }

  // We're no longer distributed if we were before
  _is_serial = true;

  // Correct our caches
  _n_nodes = 0;
  _n_elem = 0;
  _max_node_id = 0;
  _max_elem_id = 0;
  _next_free_local_node_id = libMesh::processor_id();
  _next_free_local_elem_id = libMesh::processor_id();
  _next_free_unpartitioned_node_id = libMesh::n_processors();
  _next_free_unpartitioned_elem_id = libMesh::n_processors();
}



#ifdef DEBUG
template <typename T>
void ParallelMesh::libmesh_assert_valid_parallel_object_ids 
  (const mapvector<T*> &objects) const
{
  // This function must be run on all processors at once
  parallel_only();

  unsigned int pmax_elem_id = this->parallel_max_elem_id();

  for (unsigned int i=0; i != pmax_elem_id; ++i)
    {
      T* obj = objects[i]; // Returns NULL if there's no map entry

      unsigned int dofid = obj && obj->valid_id() ?
        obj->id() : DofObject::invalid_id;
      // Local lookups by id should return the requested object
      libmesh_assert(!obj || obj->id() == i);

      unsigned int min_dofid = dofid;
      Parallel::min(min_dofid);
      // All processors with an object should agree on id
      libmesh_assert (!obj || dofid == min_dofid);

      unsigned int procid = obj && obj->valid_processor_id() ?
        obj->processor_id() : DofObject::invalid_processor_id;

      unsigned int min_procid = procid;
      Parallel::min(min_procid);
      // All processors with an object should agree on processor id
      libmesh_assert (!obj || procid == min_procid);
      // I should own any object another processor thinks I do
      libmesh_assert (min_procid != libMesh::processor_id() || obj);
    }
}



void ParallelMesh::libmesh_assert_valid_parallel_ids () const
{
  this->libmesh_assert_valid_parallel_object_ids (this->_elements);
  this->libmesh_assert_valid_parallel_object_ids (this->_nodes);
}



void ParallelMesh::libmesh_assert_valid_parallel_flags () const
{
#ifdef ENABLE_AMR
  // This function must be run on all processors at once
  parallel_only();

  unsigned int pmax_elem_id = this->parallel_max_elem_id();

  for (unsigned int i=0; i != pmax_elem_id; ++i)
    {
      Elem* elem = _elements[i]; // Returns NULL if there's no map entry

      unsigned int refinement_flag   = elem ?
        static_cast<unsigned int> (elem->refinement_flag()) : libMesh::invalid_uint;
      unsigned int p_refinement_flag = elem ?
        static_cast<unsigned int> (elem->p_refinement_flag()) : libMesh::invalid_uint;

      unsigned int min_rflag = refinement_flag;
      Parallel::min(min_rflag);
      // All processors with this element should agree on flag
      libmesh_assert (!elem || min_rflag == refinement_flag);

      unsigned int min_pflag = p_refinement_flag;
      // All processors with this element should agree on flag
      libmesh_assert (!elem || min_pflag == p_refinement_flag);
    }
#endif // ENABLE_AMR
}
#endif // DEBUG



template <typename T>
unsigned int ParallelMesh::renumber_dof_objects (mapvector<T*> &objects)
{
  // This function must be run on all processors at once
  parallel_only();

  typedef typename mapvector<T*>::veclike_iterator object_iterator;

  // In parallel we may not know what objects other processors have.
  // Start by figuring out how many
  unsigned int unpartitioned_objects = 0;

  std::vector<unsigned int>
    ghost_objects_from_proc(libMesh::n_processors(), 0);

  object_iterator it  = objects.begin();
  object_iterator end = objects.end();

  for (; it != end;)
    {
      T *obj = *it;

      // Remove any NULL container entries while we're here,
      // being careful not to invalidate our iterator
      if (!*it)
        objects.erase(it++);
      else
        {
          unsigned int obj_procid = obj->processor_id();
          if (obj_procid == DofObject::invalid_processor_id)
            unpartitioned_objects++;
          else
            ghost_objects_from_proc[obj_procid]++;
          ++it;
        }
    }

  std::vector<unsigned int> objects_on_proc(libMesh::n_processors(), 0);
  Parallel::allgather(ghost_objects_from_proc[libMesh::processor_id()],
                      objects_on_proc);

#ifndef NDEBUG
  unsigned int global_unpartitioned_objects = unpartitioned_objects;
  Parallel::max(global_unpartitioned_objects);
  libmesh_assert(global_unpartitioned_objects == unpartitioned_objects);
  for (unsigned int p=0; p != libMesh::n_processors(); ++p)
    libmesh_assert(ghost_objects_from_proc[p] <= objects_on_proc[p]);
#endif

  // We'll renumber objects in blocks by processor id
  std::vector<unsigned int> first_object_on_proc(libMesh::n_processors());
  for (unsigned int i=1; i != libMesh::n_processors(); ++i)
    first_object_on_proc[i] = first_object_on_proc[i-1] +
                              objects_on_proc[i-1];
  unsigned int next_id = first_object_on_proc[libMesh::processor_id()];
  unsigned int first_free_id =
    first_object_on_proc[libMesh::n_processors()-1] +
    objects_on_proc[libMesh::n_processors()-1] + 
    unpartitioned_objects;

  // First set new local object ids and build request sets 
  // for non-local object ids
  
  // Request sets to send to each processor
  std::vector<std::vector<unsigned int> >
    requested_ids(libMesh::n_processors());

  // We know how many objects live on each processor, so reseve() space for
  // each.
  for (unsigned int p=0; p != libMesh::n_processors(); ++p)
    if (p != libMesh::processor_id())
      requested_ids[p].reserve(ghost_objects_from_proc[p]);

  end = objects.end();
  for (it = objects.begin(); it != end; ++it)
    {
      T *obj = *it;
      if (obj->processor_id() == libMesh::processor_id())
        obj->set_id(next_id++);
      else if (obj->processor_id() != DofObject::invalid_processor_id)
        requested_ids[obj->processor_id()].push_back(obj->id());
    }

  // Next set ghost object ids from other processors
  if (libMesh::n_processors() > 1)
    {
      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();
          std::vector<unsigned int> request_to_fill;
          Parallel::send_receive(procup, requested_ids[procup],
                                 procdown, request_to_fill);

          // Fill those requests
          std::vector<unsigned int> new_ids(request_to_fill.size());
          for (unsigned int i=0; i != request_to_fill.size(); ++i)
            {
              T *obj = objects[request_to_fill[i]];
              libmesh_assert(obj);
              libmesh_assert(obj->processor_id() == libMesh::processor_id());
              new_ids[i] = obj->id();
              libmesh_assert(new_ids[i] >=
                     first_object_on_proc[libMesh::processor_id()]);
              libmesh_assert(new_ids[i] <
                     first_object_on_proc[libMesh::processor_id()] +
                     objects_on_proc[libMesh::processor_id()]);
            }

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

          // And copy the id changes we've now been informed of
          for (unsigned int i=0; i != filled_request.size(); ++i)
            {
              T *obj = objects[requested_ids[procup][i]];
              libmesh_assert (obj);
              libmesh_assert (obj->processor_id() == procup);
              libmesh_assert(filled_request[i] >=
                     first_object_on_proc[procup]);
              libmesh_assert(filled_request[i] <
                     first_object_on_proc[procup] +
                     objects_on_proc[procup]);
              obj->set_id(filled_request[i]);
            }
        }
    }

  // Next set unpartitioned object ids
  next_id = 0;
  for (unsigned int i=0; i != libMesh::n_processors(); ++i)
    next_id += objects_on_proc[i];
  for (it = objects.begin(); it != end; ++it)
    {
      T *obj = *it;
      if (obj->processor_id() == DofObject::invalid_processor_id)
        obj->set_id(next_id++);
    }

  // Finally shuffle around objects so that container indices
  // match ids
  end = objects.end();
  for (it = objects.begin(); it != end;)
    {
      T *obj = *it;
      if (obj) // don't try shuffling already-NULL entries
        {
          T *next = objects[obj->id()];
          // If we have to move this object
          if (next != obj)
            {
              // NULL out its original position for now
              // (our shuffling may put another object there shortly)
              *it = NULL;

              // There may already be another object with this id that
              // needs to be moved itself
              while (next)
                {
                  // We shouldn't be trying to give two objects the
                  // same id
                  libmesh_assert (next->id() != obj->id());
                  objects[obj->id()] = obj;
                  obj = next;
                  next = objects[obj->id()];
                }
              objects[obj->id()] = obj;
            }
        }
      // Remove any container entries that were left as NULL,
      // being careful not to invalidate our iterator
      if (!*it)
        objects.erase(it++);
      else
        ++it;
    }

  return first_free_id;
}


void ParallelMesh::renumber_nodes_and_elements ()
{
  START_LOG("renumber_nodes_and_elements()", "ParallelMesh");

#ifdef DEBUG
// Make sure our ids and flags are consistent
  this->libmesh_assert_valid_parallel_ids();
  this->libmesh_assert_valid_parallel_flags();
#endif

  std::set<unsigned int> used_nodes;

  // flag the nodes we need
  {      
    element_iterator  it = elements_begin();
    element_iterator end = elements_end();

    for (; it != end; ++it)
      {
        Elem *elem = *it;

        for (unsigned int n=0; n != elem->n_nodes(); ++n)
          used_nodes.insert(elem->node(n));
      }
  }

  // Nodes not connected to any local elements are deleted
  {
    node_iterator_imp  it = _nodes.begin();
    node_iterator_imp end = _nodes.end();

    for (; it != end;)
      {
	Node *node = *it;
        if (!used_nodes.count(node->id()))
          {
	    // remove any boundary information associated with
	    // this node
	    this->boundary_info->remove (node);

	    // delete the node
	    delete node;

            _nodes.erase(it++);
	  }
        else
          ++it;
      }
  }

  // Finally renumber all the elements
  _n_elem = this->renumber_dof_objects (this->_elements);
  _max_elem_id = _n_elem;
  _next_free_local_elem_id = _n_elem;

  // and all the remaining nodes
  _n_nodes = this->renumber_dof_objects (this->_nodes);
  _max_node_id = _n_nodes;
  _next_free_local_node_id = _n_nodes;

  // And figure out what IDs we should use when adding new nodes and
  // new elements
  unsigned int cycle = libMesh::n_processors()+1;
  unsigned int offset = _next_free_local_elem_id % cycle;
  if (offset)
    _next_free_local_elem_id += cycle - offset;
  _next_free_unpartitioned_elem_id = _next_free_local_elem_id +
                                     libMesh::n_processors();
  _next_free_local_elem_id += libMesh::processor_id();

  offset = _next_free_local_node_id % cycle;
  if (offset)
    _next_free_local_node_id += cycle - offset;
  _next_free_unpartitioned_node_id = _next_free_local_node_id +
                                     libMesh::n_processors();
  _next_free_local_node_id += libMesh::processor_id();

// Make sure our caches are up to date and our
// DofObjects are well packed
#ifdef DEBUG
  libmesh_assert(this->n_nodes() == this->parallel_n_nodes());
  libmesh_assert(this->n_elem() == this->parallel_n_elem());
  libmesh_assert(this->max_node_id() == this->parallel_max_node_id());
  libmesh_assert(this->max_elem_id() == this->parallel_max_elem_id());
  libmesh_assert(this->n_nodes() == this->max_node_id());
  libmesh_assert(this->n_elem() == this->max_elem_id());

// Make sure our ids and flags are consistent
  this->libmesh_assert_valid_parallel_ids();
  this->libmesh_assert_valid_parallel_flags();

// And make sure we've made our numbering monotonic
  MeshTools::libmesh_assert_valid_elem_ids(*this);
#endif

  STOP_LOG("renumber_nodes_and_elements()", "ParallelMesh");
}



void ParallelMesh::delete_remote_elements()
{
#ifdef DEBUG
// Make sure our neighbor links are all fine
  MeshTools::libmesh_assert_valid_neighbors(*this);
#endif

  _is_serial = false;
  MeshCommunication().delete_remote_elements(*this);

#ifdef DEBUG
// Make sure our caches are up to date and our
// DofObjects are well packed
  libmesh_assert(this->n_nodes() == this->parallel_n_nodes());
  libmesh_assert(this->n_elem() == this->parallel_n_elem());
  libmesh_assert(this->max_node_id() == this->parallel_max_node_id());
  libmesh_assert(this->max_elem_id() == this->parallel_max_elem_id());
  libmesh_assert(this->n_nodes() == this->max_node_id());
  libmesh_assert(this->n_elem() == this->max_elem_id());

// Make sure our neighbor links are all fine
  MeshTools::libmesh_assert_valid_neighbors(*this);

// Make sure our ids and flags are consistent
  this->libmesh_assert_valid_parallel_ids();
  this->libmesh_assert_valid_parallel_flags();

// And make sure our numbering is still monotonic
  MeshTools::libmesh_assert_valid_elem_ids(*this);
#endif
}



void ParallelMesh::allgather()
{
  if (_is_serial)
    return;
  _is_serial = true;
  MeshCommunication().allgather(*this);

// Make sure our caches are up to date and our
// DofObjects are well packed
#ifdef DEBUG
  libmesh_assert(this->n_nodes() == this->parallel_n_nodes());
  libmesh_assert(this->n_elem() == this->parallel_n_elem());
  libmesh_assert(this->max_node_id() == this->parallel_max_node_id());
  libmesh_assert(this->max_elem_id() == this->parallel_max_elem_id());
  libmesh_assert(this->n_nodes() == this->max_node_id());
  libmesh_assert(this->n_elem() == this->max_elem_id());

// Make sure our neighbor links are all fine
  MeshTools::libmesh_assert_valid_neighbors(*this);

// Make sure our ids and flags are consistent
  this->libmesh_assert_valid_parallel_ids();
  this->libmesh_assert_valid_parallel_flags();
#endif
}