mesh_refinement_flagging.c

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

// $Id: mesh_refinement_flagging.C 2806 2008-04-24 20:32:21Z roystgnr $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
  
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



// Local includes
#include "libmesh_config.h"

// only compile these functions if the user requests AMR support
#ifdef ENABLE_AMR

// C++ includes
#include <algorithm> // for std::sort

// Local includes
#include "elem.h"
#include "error_vector.h"
#include "mesh_refinement.h"
#include "mesh_base.h"
#include "parallel.h"



//-----------------------------------------------------------------
// Mesh refinement methods
void MeshRefinement::flag_elements_by_error_fraction (const ErrorVector& error_per_cell,
						      const Real refine_frac,
						      const Real coarsen_frac,
						      const unsigned int max_l)
{
  // 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.);

  // Clean up the refinement flags.  These could be left
  // over from previous refinement steps.
  this->clean_refinement_flags();
  
  // We're getting the minimum and maximum error values
  // for the ACTIVE elements
  Real error_min = 1.e30;
  Real error_max = 0.;

  // And, if necessary, for their parents
  Real parent_error_min = 1.e30;
  Real parent_error_max = 0.;

  // Prepare another error vector if we need to sum parent errors
  ErrorVector error_per_parent;
  if (_coarsen_by_parents) 
  {
    create_parent_error_vector(error_per_cell,
			       error_per_parent,
			       parent_error_min,
			       parent_error_max);
  }

  // We need to loop over all active elements to find the minimum 
  MeshBase::element_iterator       el_it  =
    _mesh.active_local_elements_begin();
  const MeshBase::element_iterator el_end =
    _mesh.active_local_elements_end(); 

  for (; el_it != el_end; ++el_it)
  {
    const unsigned int id  = (*el_it)->id();
    libmesh_assert (id < error_per_cell.size());

    error_max = std::max (error_max, error_per_cell[id]);
    error_min = std::min (error_min, error_per_cell[id]);
  }
  Parallel::max(error_max);
  Parallel::min(error_min);
  
  // Compute the cutoff values for coarsening and refinement
  const Real error_delta = (error_max - error_min);
  const Real parent_error_delta = parent_error_max - parent_error_min;

  const Real refine_cutoff  = (1.- _refine_fraction)*error_max;
  const Real coarsen_cutoff = _coarsen_fraction*error_delta + error_min;
  const Real parent_cutoff = _coarsen_fraction*parent_error_delta + error_min;

//   // Print information about the error
//   std::cout << " Error Information:"                     << std::endl
// 	    << " ------------------"                     << std::endl
// 	    << "   min:              " << error_min      << std::endl
// 	    << "   max:              " << error_max      << std::endl
// 	    << "   delta:            " << error_delta    << std::endl
// 	    << "     refine_cutoff:  " << refine_cutoff  << std::endl
// 	    << "     coarsen_cutoff: " << coarsen_cutoff << std::endl;
  
  

  // Loop over the elements and flag them for coarsening or
  // refinement based on the element error

  MeshBase::element_iterator       e_it  =
    _mesh.active_elements_begin();
  const MeshBase::element_iterator e_end =
    _mesh.active_elements_end(); 
  for (; e_it != e_end; ++e_it)
  {
    Elem* elem             = *e_it;
    const unsigned int id  = elem->id();

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

    if (_coarsen_by_parents)
    {
      Elem* parent           = elem->parent();
      if (parent)
      {
	const unsigned int parentid  = parent->id();
	if (error_per_parent[parentid] >= 0. &&
	    error_per_parent[parentid] <= parent_cutoff)
	  elem->set_refinement_flag(Elem::COARSEN);
      }
    }
    // Flag the element for coarsening if its error
    // is <= coarsen_fraction*delta + error_min
    else if (elem_error <= coarsen_cutoff)
    {
      elem->set_refinement_flag(Elem::COARSEN);
    }
      
    // Flag the element for refinement if its error
    // is >= refinement_cutoff.
    if (elem_error >= refine_cutoff)
      if (elem->level() < _max_h_level)
	elem->set_refinement_flag(Elem::REFINE);
  }
}



void MeshRefinement::flag_elements_by_error_tolerance (const ErrorVector& error_per_cell_in)
{
  // Check for valid fractions..
  // The fraction values must be in [0,1]
  libmesh_assert (_coarsen_threshold  >= 0. && _refine_fraction  <= 1.);
  libmesh_assert (_refine_fraction  >= 0. && _refine_fraction  <= 1.);
  libmesh_assert (_coarsen_fraction >= 0. && _coarsen_fraction <= 1.);

  // How much error per cell will we tolerate?
  const Real local_refinement_tolerance =
    _absolute_global_tolerance / std::sqrt(static_cast<Real>(_mesh.n_active_elem()));
  const Real local_coarsening_tolerance =
    local_refinement_tolerance * _coarsen_threshold;

  // Prepare another error vector if we need to sum parent errors
  ErrorVector error_per_parent;
  if (_coarsen_by_parents) 
  {
    Real parent_error_min, parent_error_max;

    create_parent_error_vector(error_per_cell_in,
			       error_per_parent,
			       parent_error_min,
			       parent_error_max);
  }

  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;
    Elem* parent = elem->parent();
    const unsigned int elem_number = elem->id();
    const float        elem_error  = error_per_cell_in[elem_number];

    if (elem_error > local_refinement_tolerance &&
	elem->level() < _max_h_level)
      elem->set_refinement_flag(Elem::REFINE);

    if (!_coarsen_by_parents && elem_error <
	local_coarsening_tolerance)
      elem->set_refinement_flag(Elem::COARSEN);

    if (_coarsen_by_parents && parent)
    {
      float parent_error = error_per_parent[parent->id()];
      if (parent_error >= 0.)
      {
	const Real parent_coarsening_tolerance =
	  std::sqrt(parent->n_children() *
		    local_coarsening_tolerance *
		    local_coarsening_tolerance);
	if (parent_error < parent_coarsening_tolerance)
	  elem->set_refinement_flag(Elem::COARSEN);
      }
    }
  }
}



bool MeshRefinement::flag_elements_by_nelem_target (const ErrorVector& error_per_cell)
{
  // 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.);

  // This function is currently only coded to work when coarsening by
  // parents - it's too hard to guess how many coarsenings will be
  // performed otherwise.
  libmesh_assert (_coarsen_by_parents);
  
  // The number of active elements in the mesh - hopefully less than
  // 2 billion on 32 bit machines
  const unsigned int n_active_elem  = _mesh.n_active_elem();

  // The maximum number of active elements to flag for coarsening
  const unsigned int max_elem_coarsen =
    static_cast<unsigned int>(_coarsen_fraction * n_active_elem) + 1;

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

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

  // The target number of elements to add or remove
  const int n_elem_new = _nelem_target - n_active_elem;

  // Create an vector with active element errors and ids,
  // sorted by highest errors first
  std::vector<std::pair<float, unsigned int> > sorted_error;

  sorted_error.reserve (n_active_elem);

  // 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_nelem_target does not yet "
                  << "work on a parallel mesh." << std::endl;
      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)
    {
      unsigned int eid = (*elem_it)->id();
      libmesh_assert(eid < error_per_cell.size());
      sorted_error.push_back
        (std::make_pair(error_per_cell[eid], eid));
    }

  // Default sort works since pairs are sorted lexicographically
  std::sort (sorted_error.begin(), sorted_error.end());
  std::reverse (sorted_error.begin(), sorted_error.end());

  // Create a sorted error vector with coarsenable parent elements
  // only, sorted by lowest errors first
  ErrorVector error_per_parent;
  std::vector<std::pair<float, unsigned int> > sorted_parent_error;
  Real parent_error_min, parent_error_max;

  create_parent_error_vector(error_per_cell,
                             error_per_parent,
                             parent_error_min,
                             parent_error_max);

  // create_parent_error_vector sets values for non-parents and
  // non-coarsenable parents to -1.  Get rid of them.
  for (unsigned int i=0; i != error_per_parent.size(); ++i)
    if (error_per_parent[i] != -1)
      sorted_parent_error.push_back(std::make_pair(error_per_parent[i], i));

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

  // Keep track of how many elements we plan to coarsen & refine
  unsigned int coarsen_count = 0;
  unsigned int refine_count = 0;

  const unsigned int dim = _mesh.mesh_dimension();
  unsigned int twotodim = 1;
  for (unsigned int i=0; i!=dim; ++i)
    twotodim *= 2;

  // First, let's try to get our element count to target_nelem
  if (n_elem_new >= 0)
  {
    // Every element refinement creates at least 
    // 2^dim-1 new elements
    refine_count =
      std::min(static_cast<unsigned int>(n_elem_new / (twotodim-1)),
	       max_elem_refine);
  }
  else
  {
    // Every successful element coarsening is likely to destroy
    // 2^dim-1 net elements.
    coarsen_count =
      std::min(static_cast<unsigned int>(-n_elem_new / (twotodim-1)),
	       max_elem_coarsen);
  }

  // Next, let's see if we can trade any refinement for coarsening