mesh_refinement.c

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

// $Id: mesh_refinement.C 2840 2008-05-12 12:56:32Z 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


// C++ includes
#include <cmath> // for isnan(), when it's defined
#include <limits>

// Local includes
#include "libmesh_config.h"

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

#include "boundary_info.h"
#include "elem.h"
#include "error_vector.h"
#include "libmesh_logging.h"
#include "mesh_base.h"
#include "mesh_communication.h"
#include "mesh_refinement.h"
#include "parallel.h"
#include "parallel_ghost_sync.h"
#include "remote_elem.h"

#ifdef DEBUG
// Some extra validation for ParallelMesh
#include "mesh_tools.h"
#include "parallel_mesh.h"
#endif // DEBUG


//-----------------------------------------------------------------
// Mesh refinement methods
MeshRefinement::MeshRefinement (MeshBase& m) :
  _mesh(m),
  _use_member_parameters(false),
  _coarsen_by_parents(false),
  _refine_fraction(0.3),
  _coarsen_fraction(0.0),
  _max_h_level(libMesh::invalid_uint),
  _coarsen_threshold(10),
  _nelem_target(0),
  _absolute_global_tolerance(0.0),
  _face_level_mismatch_limit(1),
  _edge_level_mismatch_limit(0),
  _node_level_mismatch_limit(0)
{
}



MeshRefinement::~MeshRefinement ()
{
  this->clear();  
}



void MeshRefinement::clear ()
{
  _new_nodes_map.clear();
}



Node* MeshRefinement::add_point (const Point& p,
                                 const unsigned int processor_id,
                                 const Real tol)
{
  START_LOG("add_point()", "MeshRefinement");

  // Return the node if it already exists
  Node *node = _new_nodes_map.find(p, tol);
  if (node)
    {
      STOP_LOG("add_point()", "MeshRefinement");
      return node;
    }
      
  // Add the node, with a default id and the requested
  // processor_id
  node = _mesh.add_point (p, DofObject::invalid_id, processor_id);

  libmesh_assert (node != NULL);

  // Add the node to the map.
  _new_nodes_map.insert(*node);

  // Return the address of the new node
  STOP_LOG("add_point()", "MeshRefinement");
  return node;
}



Elem* MeshRefinement::add_elem (Elem* elem)
{
  libmesh_assert (elem != NULL);

  
//   // If the unused_elements has any iterators from
//   // old elements, take the first one
//   if (!_unused_elements.empty())
//     {
//       std::vector<Elem*>::iterator it = _unused_elements.front();

//       *it = elem;

//       _unused_elements.pop_front();
//     }

//   // Otherwise, use the conventional add method
//   else
//     {
//       _mesh.add_elem (elem);
//     }

  // The _unused_elements optimization has been turned off.
  _mesh.add_elem (elem);
  
  return elem;
}



void MeshRefinement::create_parent_error_vector
  (const ErrorVector& error_per_cell,
   ErrorVector& error_per_parent,
   Real& parent_error_min,
   Real& parent_error_max)
{
  // This function must be run on all processors at once
  parallel_only();

  // Make sure the input error vector is valid
#ifdef DEBUG
  for (unsigned int i=0; i != error_per_cell.size(); ++i)
    {
      libmesh_assert(error_per_cell[i] >= 0);
  // isnan() isn't standard C++ yet
  #ifdef isnan
      libmesh_assert(!isnan(error_per_cell[i]));
  #endif
    }
#endif // #ifdef DEBUG

  // error values on uncoarsenable elements will be left at -1
  error_per_parent.clear();
  error_per_parent.resize(error_per_cell.size(), 0.0);

  {
  // Find which elements are uncoarsenable
  MeshBase::element_iterator       elem_it  = _mesh.active_local_elements_begin();
  const MeshBase::element_iterator elem_end = _mesh.active_local_elements_end();
  for (; elem_it != elem_end; ++elem_it)
    {
      Elem* elem   = *elem_it;
      Elem* parent = elem->parent();

      // Active elements are uncoarsenable
      error_per_parent[elem->id()] = -1.0;

      // Grandparents and up are uncoarsenable
      while (parent)
        {
          parent = parent->parent();
          if (parent)
            {
              const unsigned int parentid  = parent->id();
              libmesh_assert (parentid < error_per_parent.size());
              error_per_parent[parentid] = -1.0;
            }
        }
    }

  // Sync between processors.
  // Use a reference to std::vector to avoid confusing
  // Parallel::min
  std::vector<ErrorVectorReal> &epp = error_per_parent;
  Parallel::min(epp);
  }

  // The parent's error is defined as the square root of the
  // sum of the children's errors squared, so errors that are
  // Hilbert norms remain Hilbert norms.
  //
  // Because the children may be on different processors, we
  // calculate local contributions to the parents' errors squared
  // first, then sum across processors and take the square roots
  // second.
  {
  MeshBase::element_iterator       elem_it  = _mesh.active_local_elements_begin();
  const MeshBase::element_iterator elem_end = _mesh.active_local_elements_end();

  for (; elem_it != elem_end; ++elem_it)
    {
      Elem* elem   = *elem_it;
      Elem* parent = elem->parent();

      // Calculate each contribution to parent cells
      if (parent)
        {
          const unsigned int parentid  = parent->id();
          libmesh_assert (parentid < error_per_parent.size());

	  // If the parent has grandchildren we won't be able to
	  // coarsen it, so forget it.  Otherwise, add this child's
	  // contribution to the sum of the squared child errors
	  if (error_per_parent[parentid] != -1.0)
            error_per_parent[parentid] += (error_per_cell[elem->id()] *
                                           error_per_cell[elem->id()]);
        }
    }
  }

  // Sum the vector across all processors
  Parallel::sum(static_cast<std::vector<ErrorVectorReal>&>(error_per_parent));

  // Calculate the min and max as we loop
  parent_error_min = std::numeric_limits<double>::max();
  parent_error_max = 0.;

  for (unsigned int i = 0; i != error_per_parent.size(); ++i)
    {
      // If this element isn't a coarsenable parent with error, we
      // have nothing to do.  Just flag it as -1 and move on
      // Note that Parallel::sum might have left uncoarsenable
      // elements with error_per_parent=-n_proc, so reset it to
      // error_per_parent=-1
      if (error_per_parent[i] < 0.)
        {
          error_per_parent[i] = -1.;
          continue;
        }

      // The error estimator might have already given us an
      // estimate on the coarsenable parent elements; if so then
      // we want to retain that estimate
      if (error_per_cell[i])
        {
          error_per_parent[i] = error_per_cell[i];
          continue;
        }
      // if not, then e_parent = sqrt(sum(e_child^2))
      else
        error_per_parent[i] = std::sqrt(error_per_parent[i]);

      parent_error_min = std::min (parent_error_min,
                                   error_per_parent[i]);
      parent_error_max = std::max (parent_error_max,
                                   error_per_parent[i]);
    }
}



void MeshRefinement::update_nodes_map ()
{
  this->_new_nodes_map.init(_mesh);
} 



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

  MeshBase::element_iterator       elem_it  = _mesh.active_local_elements_begin();
  const MeshBase::element_iterator elem_end = _mesh.active_local_elements_end();

  bool failure = false;

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

      for (unsigned int n=0; n<elem->n_neighbors(); n++)
        {
          Elem *neighbor = elem->neighbor(n);

          if (!neighbor || !neighbor->active() ||
              neighbor == remote_elem)
            continue;

          if ((neighbor->level() + 1 < elem->level()) ||
              (neighbor->p_level() + 1 < elem->p_level()) ||
              (neighbor->p_level() > elem->p_level() + 1))
            {
              failure = true;
              break;
            }
        }
    }

  // If any processor failed, we failed globally
  Parallel::max(failure);

  if (failure)
    {
      // We didn't pass the level one test, so libmesh_assert that 
      // we're allowed not to
      libmesh_assert(!libmesh_assert_pass);
      return false;
    }
  return true;
}



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

  bool found_flag = false;

  // Search for local flags
  MeshBase::element_iterator       elem_it  = _mesh.active_local_elements_begin();
  const MeshBase::element_iterator elem_end = _mesh.active_local_elements_end();

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

      if (elem->refinement_flag() == Elem::REFINE ||
          elem->refinement_flag() == Elem::COARSEN ||
          elem->p_refinement_flag() == Elem::REFINE ||
          elem->p_refinement_flag() == Elem::COARSEN)
        {
          found_flag = true;
          break;
        }
    }

  // If we found a flag on any processor, it counts
  Parallel::max(found_flag);

  if (found_flag)
    {
      // We didn't pass the "elements are unflagged" test,
      // so libmesh_assert that we're allowed not to
      libmesh_assert(!libmesh_assert_pass);
      return false;
    }
  return true;
}



bool MeshRefinement::refine_and_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);
    }