mesh_refinement_flagging.c

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

  while (coarsen_count < max_elem_coarsen &&
         refine_count < max_elem_refine &&
         coarsen_count < sorted_parent_error.size() &&
         refine_count < sorted_error.size() &&
         sorted_error[refine_count].first > 
	 sorted_parent_error[coarsen_count].first * _coarsen_threshold)
  {
    coarsen_count++;
    refine_count++;
  }

  if (refine_count > max_elem_refine)
    refine_count = max_elem_refine;
  unsigned int successful_refine_count = 0;
  for (unsigned int i=0; i != sorted_error.size(); ++i)
    {
      if (successful_refine_count >= refine_count)
        break;

      unsigned int eid = sorted_error[i].second;
      Elem *elem = _mesh.elem(eid);
      if (elem->level() < _max_h_level)
        {
	  elem->set_refinement_flag(Elem::REFINE);
	  successful_refine_count++;
        }
    }

  // If we couldn't refine enough elements, don't coarsen too many
  // either
  if (coarsen_count < (refine_count - successful_refine_count))
    coarsen_count = 0;
  else
    coarsen_count -= (refine_count - successful_refine_count);

  if (coarsen_count > max_elem_coarsen)
    coarsen_count = max_elem_coarsen;
  unsigned int successful_coarsen_count = 0;
  for (unsigned int i=0; i != sorted_parent_error.size(); ++i)
    {
      if (successful_coarsen_count >= coarsen_count * twotodim)
        break;

      unsigned int parent_id = sorted_parent_error[i].second;
      Elem *parent = _mesh.elem(parent_id);
      libmesh_assert(parent->has_children());
      for (unsigned int c=0; c != parent->n_children(); ++c)
        {
          Elem *elem = parent->child(c);
          if (elem->active())
            {
              elem->set_refinement_flag(Elem::COARSEN);
              successful_coarsen_count++;
            }
        }
    }

  // Return true if we've done all the AMR/C we can
  if (!successful_coarsen_count && 
      !successful_refine_count)
    return true;
  // And false if there may still be more to do.
  return false;
}



void MeshRefinement::flag_elements_by_elem_fraction (const ErrorVector& error_per_cell,
						     const Real refine_frac,
						     const Real coarsen_frac,
						     const unsigned int max_l)
{
  // FIXME - this won't work on a non-serialized mesh yet
  if (!_mesh.is_serial())
    {
      if (libMesh::processor_id() == 0)
        std::cerr << "flag_elements_by_elem_fraction does not yet "
                  << "work on a parallel mesh." << std::endl;
      error();
    }

  // The function arguments are currently just there for
  // backwards_compatibility
  if (!_use_member_parameters)
  {
    // If the user used non-default parameters, lets warn
    // that they're deprecated
    if (refine_frac != 0.3 ||
	coarsen_frac != 0.0 ||
	max_l != libMesh::invalid_uint)
      deprecated();

    _refine_fraction = refine_frac;
    _coarsen_fraction = coarsen_frac;
    _max_h_level = max_l;
  }

  // Check for valid fractions..
  // The fraction values must be in [0,1]
  libmesh_assert (_refine_fraction  >= 0. && _refine_fraction  <= 1.);
  libmesh_assert (_coarsen_fraction >= 0. && _coarsen_fraction <= 1.);

  // The number of active elements in the mesh
  const unsigned int n_active_elem  = _mesh.n_elem();

  // The number of elements to flag for coarsening
  const unsigned int n_elem_coarsen =
    static_cast<unsigned int>(_coarsen_fraction * n_active_elem);

  // The number of elements to flag for refinement
  const unsigned int n_elem_refine =
    static_cast<unsigned int>(_refine_fraction  * n_active_elem);


  
  // Clean up the refinement flags.  These could be left
  // over from previous refinement steps.
  this->clean_refinement_flags();


  // This vector stores the error and element number for all the
  // active elements.  It will be sorted and the top & bottom
  // elements will then be flagged for coarsening & refinement
  std::vector<float> sorted_error;

  sorted_error.reserve (n_active_elem);

  // Loop over the active elements and create the entry
  // in the sorted_error vector
  MeshBase::element_iterator       elem_it  = _mesh.active_elements_begin();
  const MeshBase::element_iterator elem_end = _mesh.active_elements_end(); 

  for (; elem_it != elem_end; ++elem_it)
    sorted_error.push_back (error_per_cell[(*elem_it)->id()]);

  // Now sort the sorted_error vector
  std::sort (sorted_error.begin(), sorted_error.end());
  
  // If we're coarsening by parents:
  // Create a sorted error vector with coarsenable parent elements
  // only, sorted by lowest errors first
  ErrorVector error_per_parent, sorted_parent_error;
  if (_coarsen_by_parents)
  {
    Real parent_error_min, parent_error_max;

    create_parent_error_vector(error_per_cell,
			       error_per_parent,
			       parent_error_min,
			       parent_error_max);

    sorted_parent_error = error_per_parent;
    std::sort (sorted_parent_error.begin(), sorted_parent_error.end());

    // All the other error values will be 0., so get rid of them.
    sorted_parent_error.erase (std::remove(sorted_parent_error.begin(),
					   sorted_parent_error.end(), 0.),
			       sorted_parent_error.end());
  }
  

  float top_error= 0., bottom_error = 0.;

  // Get the maximum error value corresponding to the
  // bottom n_elem_coarsen elements
  if (_coarsen_by_parents && n_elem_coarsen)
    {
      const unsigned int dim = _mesh.mesh_dimension();
      unsigned int twotodim = 1;
      for (unsigned int i=0; i!=dim; ++i)
        twotodim *= 2;

      unsigned int n_parent_coarsen = n_elem_coarsen / (twotodim - 1);

      if (n_parent_coarsen)
	bottom_error = sorted_parent_error[n_parent_coarsen - 1];
    }
  else if (n_elem_coarsen)
    {
      bottom_error = sorted_error[n_elem_coarsen - 1];
    }

  if (n_elem_refine)
    top_error = sorted_error[sorted_error.size() - n_elem_refine + 1];

  // Finally, let's do the element flagging
  elem_it  = _mesh.active_elements_begin();
  for (; elem_it != elem_end; ++elem_it)
    {
      Elem* elem = *elem_it;
      Elem* parent = elem->parent();

      if (_coarsen_by_parents && parent && n_elem_coarsen &&
          error_per_parent[parent->id()] <= bottom_error)
        elem->set_refinement_flag(Elem::COARSEN);

      if (!_coarsen_by_parents && n_elem_coarsen &&
          error_per_cell[elem->id()] <= bottom_error)
        elem->set_refinement_flag(Elem::COARSEN);

      if (n_elem_refine &&
          elem->level() < _max_h_level &&
          error_per_cell[elem->id()] >= top_error)
        elem->set_refinement_flag(Elem::REFINE);
    }
}



void MeshRefinement::flag_elements_by_mean_stddev (const ErrorVector& error_per_cell,
						   const Real refine_frac,
						   const Real coarsen_frac,
						   const unsigned int max_l)
{
  // FIXME - this won't work on a non-serialized mesh yet
  if (!_mesh.is_serial())
    {
      if (libMesh::processor_id() == 0)
        std::cerr << "flag_elements_by_mean_stddev does not yet "
                  << "work on a parallel mesh." << std::endl;
      error();
    }

  // The function arguments are currently just there for
  // backwards_compatibility
  if (!_use_member_parameters)
    {
      // If the user used non-default parameters, lets warn
      // that they're deprecated
      if (refine_frac != 0.3 ||
          coarsen_frac != 0.0 ||
          max_l != libMesh::invalid_uint)
        deprecated();

      _refine_fraction = refine_frac;
      _coarsen_fraction = coarsen_frac;
      _max_h_level = max_l;
    }

  // Get the mean value from the error vector
  const Real mean = error_per_cell.mean();
  
  // Get the standard deviation.  This equals the
  // square-root of the variance
  const Real stddev = std::sqrt (error_per_cell.variance());
  
  // Check for valid fractions
  libmesh_assert (_refine_fraction  >= 0. && _refine_fraction  <= 1.);
  libmesh_assert (_coarsen_fraction >= 0. && _coarsen_fraction <= 1.);

  // The refine and coarsen cutoff
  const Real refine_cutoff  =  mean + _refine_fraction  * stddev;
  const Real coarsen_cutoff =  std::max(mean - _coarsen_fraction * stddev, 0.);
      
  // Loop over the elements and flag them for coarsening or
  // refinement based on the element error
  MeshBase::element_iterator       elem_it  = _mesh.active_elements_begin();
  const MeshBase::element_iterator elem_end = _mesh.active_elements_end(); 

  for (; elem_it != elem_end; ++elem_it)
    {
      Elem* elem             = *elem_it;
      const unsigned int id  = elem->id();

      libmesh_assert (id < error_per_cell.size());
      
      const float elem_error = error_per_cell[id];

      // Possibly flag the element for coarsening ...
      if (elem_error <= coarsen_cutoff)
	elem->set_refinement_flag(Elem::COARSEN);
      
      // ... or refinement
      if ((elem_error >= refine_cutoff) && (elem->level() < _max_h_level))
	elem->set_refinement_flag(Elem::REFINE);
    }
}



void MeshRefinement::switch_h_to_p_refinement ()
{
  MeshBase::element_iterator       elem_it  = _mesh.elements_begin();
  const MeshBase::element_iterator elem_end = _mesh.elements_end(); 

  for ( ; elem_it != elem_end; ++elem_it)
    {
      if ((*elem_it)->active())
        {
          (*elem_it)->set_p_refinement_flag((*elem_it)->refinement_flag());
          (*elem_it)->set_refinement_flag(Elem::DO_NOTHING);
        } 
      else
        {
          (*elem_it)->set_p_refinement_flag((*elem_it)->refinement_flag());
          (*elem_it)->set_refinement_flag(Elem::INACTIVE);
        } 
    }
}



void MeshRefinement::add_p_to_h_refinement ()
{
  MeshBase::element_iterator       elem_it  = _mesh.elements_begin();
  const MeshBase::element_iterator elem_end = _mesh.elements_end(); 

  for ( ; elem_it != elem_end; ++elem_it)
    (*elem_it)->set_p_refinement_flag((*elem_it)->refinement_flag());
}



void MeshRefinement::clean_refinement_flags ()
{
  // Possibly clean up the refinement flags from
  // a previous step
//   elem_iterator       elem_it (_mesh.elements_begin());
//   const elem_iterator elem_end(_mesh.elements_end());

  MeshBase::element_iterator       elem_it  = _mesh.elements_begin();
  const MeshBase::element_iterator elem_end = _mesh.elements_end(); 

  for ( ; elem_it != elem_end; ++elem_it)
    {
      if ((*elem_it)->active())
        {
          (*elem_it)->set_refinement_flag(Elem::DO_NOTHING);
          (*elem_it)->set_p_refinement_flag(Elem::DO_NOTHING);
        } 
      else
        {
          (*elem_it)->set_refinement_flag(Elem::INACTIVE);
          (*elem_it)->set_p_refinement_flag(Elem::INACTIVE);
        } 
    }
}

#endif