mesh_refinement.c

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

  // Parallel consistency has to come first, or coarsening
  // along processor boundaries might occasionally be falsely
  // prevented
  if (!_mesh.is_serial())
    this->make_flags_parallel_consistent();

  // Repeat until flag changes match on every processor
  do
    {
      // Repeat until coarsening & refinement flags jive
      bool satisfied = false;
      do
        {
          const bool coarsening_satisfied =
	    this->make_coarsening_compatible(maintain_level_one);

          const bool refinement_satisfied =
	    this->make_refinement_compatible(maintain_level_one);

          bool smoothing_satisfied = 
 	    !this->eliminate_unrefined_patches();

          if (_edge_level_mismatch_limit)
            smoothing_satisfied = smoothing_satisfied && 
              !this->limit_level_mismatch_at_edge (_edge_level_mismatch_limit);

          if (_node_level_mismatch_limit)
            smoothing_satisfied = smoothing_satisfied &&
              !this->limit_level_mismatch_at_node (_node_level_mismatch_limit);

          satisfied = (coarsening_satisfied &&
		       refinement_satisfied &&
		       smoothing_satisfied);
#ifdef DEBUG
          bool max_satisfied = satisfied,
               min_satisfied = satisfied;
          Parallel::max(max_satisfied);
          Parallel::min(min_satisfied);
          libmesh_assert (satisfied == max_satisfied);
          libmesh_assert (satisfied == min_satisfied);
#endif
        }
      while (!satisfied);
    }
  while (!_mesh.is_serial() && !this->make_flags_parallel_consistent());

  // First coarsen the flagged elements.
  const bool coarsening_changed_mesh =
    this->_coarsen_elements ();

  if (_maintain_level_one)
    libmesh_assert(test_level_one(true));
  libmesh_assert(this->make_coarsening_compatible(maintain_level_one));
  libmesh_assert(this->make_refinement_compatible(maintain_level_one));
// FIXME: This won't pass unless we add a redundant find_neighbors()
// call or replace find_neighbors() with on-the-fly neighbor updating
// libmesh_assert(!this->eliminate_unrefined_patches());

  // We can't contract the mesh ourselves anymore - a System might
  // need to restrict old coefficient vectors first
  // _mesh.contract();

  // Now refine the flagged elements.  This will
  // take up some space, maybe more than what was freed.
  const bool refining_changed_mesh =
    this->_refine_elements();

  if (_maintain_level_one)
    libmesh_assert(test_level_one(true));
  libmesh_assert(test_unflagged(true));
  libmesh_assert(this->make_coarsening_compatible(maintain_level_one));
  libmesh_assert(this->make_refinement_compatible(maintain_level_one));
// FIXME: This won't pass unless we add a redundant find_neighbors()
// call or replace find_neighbors() with on-the-fly neighbor updating
// libmesh_assert(!this->eliminate_unrefined_patches());

  // Finally, the new mesh needs to be prepared for use
  if (coarsening_changed_mesh || refining_changed_mesh)
    {
#ifdef DEBUG
      ParallelMesh *pmesh = dynamic_cast<ParallelMesh *>(&_mesh);
      if (pmesh)
        pmesh->libmesh_assert_valid_parallel_ids();
#endif
      _mesh.prepare_for_use ();
      
      return true;
    }

  // Otherwise there was no change in the mesh,
  // let the user know.  Also, there is no need
  // to prepare the mesh for use since it did not change.
  return false;
    
}







bool MeshRefinement::coarsen_elements (const bool maintain_level_one)
{
  // This function must be run on all processors at once
  parallel_only();

  bool _maintain_level_one = maintain_level_one;

  // If the user used non-default parameters, let's warn that they're
  // deprecated
  if (!maintain_level_one)
    {
      deprecated();
    }
  else
    _maintain_level_one = _face_level_mismatch_limit;

  // We can't yet turn a non-level-one mesh into a level-one mesh
  if (_maintain_level_one)
    libmesh_assert(test_level_one(true));

  // Possibly clean up the refinement flags from
  // a previous step
  MeshBase::element_iterator       elem_it  = _mesh.elements_begin();
  const MeshBase::element_iterator elem_end = _mesh.elements_end();

  for ( ; elem_it != elem_end; ++elem_it)
    {
      // Pointer to the element
      Elem* elem = *elem_it;

      // Set refinement flag to INACTIVE if the
      // element isn't active
      if ( !elem->active())
        {
	  elem->set_refinement_flag(Elem::INACTIVE);
	  elem->set_p_refinement_flag(Elem::INACTIVE);
        }

      // This might be left over from the last step
      if (elem->refinement_flag() == Elem::JUST_REFINED)
	elem->set_refinement_flag(Elem::DO_NOTHING);
    }

  // Parallel consistency has to come first, or coarsening
  // along processor boundaries might occasionally be falsely
  // prevented
  if (!_mesh.is_serial())
    this->make_flags_parallel_consistent();

  // Repeat until flag changes match on every processor
  do
    {
      // Repeat until the flags form a conforming mesh.
      bool satisfied = false;
      do
        {
          const bool coarsening_satisfied =
	    this->make_coarsening_compatible(maintain_level_one);

          bool smoothing_satisfied = 
 	    !this->eliminate_unrefined_patches();// &&

          if (_edge_level_mismatch_limit)
            smoothing_satisfied = smoothing_satisfied &&
              !this->limit_level_mismatch_at_edge (_edge_level_mismatch_limit);

          if (_node_level_mismatch_limit)
            smoothing_satisfied = smoothing_satisfied &&
              !this->limit_level_mismatch_at_node (_node_level_mismatch_limit);

          satisfied = (coarsening_satisfied &&
		       smoothing_satisfied);
#ifdef DEBUG
          bool max_satisfied = satisfied,
               min_satisfied = satisfied;
          Parallel::max(max_satisfied);
          Parallel::min(min_satisfied);
          libmesh_assert (satisfied == max_satisfied);
          libmesh_assert (satisfied == min_satisfied);
#endif
        }
      while (!satisfied);
    }
  while (!_mesh.is_serial() && !this->make_flags_parallel_consistent());

  // Coarsen the flagged elements.
  const bool mesh_changed = 
    this->_coarsen_elements ();

  if (_maintain_level_one)
    libmesh_assert(test_level_one(true));
  libmesh_assert(this->make_coarsening_compatible(maintain_level_one));
// FIXME: This won't pass unless we add a redundant find_neighbors()
// call or replace find_neighbors() with on-the-fly neighbor updating
// libmesh_assert(!this->eliminate_unrefined_patches());
    
  // We can't contract the mesh ourselves anymore - a System might
  // need to restrict old coefficient vectors first
  // _mesh.contract();

  // Finally, the new mesh may need to be prepared for use
  if (mesh_changed)
    _mesh.prepare_for_use ();

  return mesh_changed;
}







bool MeshRefinement::refine_elements (const bool maintain_level_one)
{
  // This function must be run on all processors at once
  parallel_only();

  bool _maintain_level_one = maintain_level_one;

  // If the user used non-default parameters, let's warn that they're
  // deprecated
  if (!maintain_level_one)
    {
      deprecated();
    }
  else
    _maintain_level_one = _face_level_mismatch_limit;

  if (_maintain_level_one)
    libmesh_assert(test_level_one(true));

  // Possibly clean up the refinement flags from
  // a previous step
  MeshBase::element_iterator       elem_it  = _mesh.elements_begin();
  const MeshBase::element_iterator elem_end = _mesh.elements_end();

  for ( ; elem_it != elem_end; ++elem_it)
    {
      // Pointer to the element
      Elem *elem = *elem_it;

      // Set refinement flag to INACTIVE if the
      // element isn't active
      if ( !elem->active())
        {
	  elem->set_refinement_flag(Elem::INACTIVE);
	  elem->set_p_refinement_flag(Elem::INACTIVE);
        }

      // This might be left over from the last step
      if (elem->refinement_flag() == Elem::JUST_REFINED)
	elem->set_refinement_flag(Elem::DO_NOTHING);
    }


  
  // Parallel consistency has to come first, or coarsening
  // along processor boundaries might occasionally be falsely
  // prevented
  if (!_mesh.is_serial())
    this->make_flags_parallel_consistent();

  // Repeat until flag changes match on every processor
  do
    {
      // Repeat until coarsening & refinement flags jive
      bool satisfied = false;
      do
        {
          const bool refinement_satisfied =
	    this->make_refinement_compatible(maintain_level_one);

          bool smoothing_satisfied = 
 	    !this->eliminate_unrefined_patches();// &&

          if (_edge_level_mismatch_limit)
            smoothing_satisfied = smoothing_satisfied &&
              !this->limit_level_mismatch_at_edge (_edge_level_mismatch_limit);

          if (_node_level_mismatch_limit)
            smoothing_satisfied = smoothing_satisfied &&
              !this->limit_level_mismatch_at_node (_node_level_mismatch_limit);

          satisfied = (refinement_satisfied &&
		       smoothing_satisfied);
#ifdef DEBUG
          bool max_satisfied = satisfied,
               min_satisfied = satisfied;
          Parallel::max(max_satisfied);
          Parallel::min(min_satisfied);
          libmesh_assert (satisfied == max_satisfied);
          libmesh_assert (satisfied == min_satisfied);
#endif
        }
      while (!satisfied);
    }
  while (!_mesh.is_serial() && !this->make_flags_parallel_consistent());
  
  // Now refine the flagged elements.  This will
  // take up some space, maybe more than what was freed.
  const bool mesh_changed = 
    this->_refine_elements();

  if (_maintain_level_one)
    libmesh_assert(test_level_one(true));
  libmesh_assert(this->make_refinement_compatible(maintain_level_one));
// FIXME: This won't pass unless we add a redundant find_neighbors()
// call or replace find_neighbors() with on-the-fly neighbor updating
// libmesh_assert(!this->eliminate_unrefined_patches());
    
  // Finally, the new mesh needs to be prepared for use
  if (mesh_changed)
    _mesh.prepare_for_use ();

  return mesh_changed;
}


// Functor for make_flags_parallel_consistent
namespace {

struct SyncRefinementFlags
{
typedef unsigned char datum;
typedef Elem::RefinementState (Elem::*get_a_flag)() const;
typedef void (Elem::*set_a_flag)(const Elem::RefinementState);

SyncRefinementFlags(MeshBase &_mesh,
                    get_a_flag _getter,
                    set_a_flag _setter) :
  mesh(_mesh), parallel_consistent(true),
  get_flag(_getter), set_flag(_setter) {}

MeshBase &mesh;
bool parallel_consistent;
get_a_flag get_flag;
set_a_flag set_flag;
// References to pointers to member functions segfault?
// get_a_flag& get_flag;
// set_a_flag& set_flag;

// Find the refinement flag on each requested element
void gather_data (const std::vector<unsigned int>& ids,
                  std::vector<datum>& flags)
{
  flags.resize(ids.size());

  for (unsigned int i=0; i != ids.size(); ++i)
    {
      // Look for this element in the mesh
      Elem *elem = mesh.elem(ids[i]);

      // We'd better find every element we're asked for
      libmesh_assert (elem);

      // Return the element's refinement flag
      flags[i] = (elem->*get_flag)();
    }
}

void act_on_data (const std::vector<unsigned int>& ids,
                  std::vector<datum>& flags)
{
  for (unsigned int i=0; i != ids.size(); ++i)
    {
      Elem *elem = mesh.elem(ids[i]);

      libmesh_assert(elem);

      datum old_flag = (elem->*get_flag)();
      datum &new_flag = flags[i];

      if (old_flag != new_flag)
        {
	  // It's possible for foreign flags to be (temporarily) more
	  // conservative than our own, such as when a refinement in
	  // one of the foreign processor's elements is mandated by a
	  // refinement in one of our neighboring elements it can see
	  // which was mandated by a refinement in one of our
	  // neighboring elements it can't see
          // libmesh_assert (!(new_flag != Elem::REFINE && 
          //                   old_flag == Elem::REFINE));
	  //
          (elem->*set_flag)
            (static_cast<Elem::RefinementState>(new_flag));
          parallel_consistent = false;
        }
    }
}
};
}



bool MeshRefinement::make_flags_parallel_consistent()
{
  // This function must be run on all processors at once
  parallel_only();

  START_LOG ("make_flags_parallel_consistent()", "MeshRefinement");

  SyncRefinementFlags hsync(_mesh, &Elem::refinement_flag,
                            &Elem::set_refinement_flag);
  Parallel::sync_dofobject_data_by_id
    (_mesh.elements_begin(), _mesh.elements_end(), hsync);

  SyncRefinementFlags psync(_mesh, &Elem::p_refinement_flag,
                            &Elem::set_p_refinement_flag);
  Parallel::sync_dofobject_data_by_id