uniform_refinement_estimator.c

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

// $Id: uniform_refinement_estimator.C 2921 2008-07-08 21:23:50Z jwpeterson $

// 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 <algorithm> // for std::fill
#include <cmath>    // for sqrt


// Local Includes
#include "dof_map.h"
#include "elem.h"
#include "equation_systems.h"
#include "error_vector.h"
#include "fe.h"
#include "fe_interface.h"
#include "libmesh_common.h"
#include "libmesh_logging.h"
#include "mesh.h"
#include "mesh_refinement.h"
#include "numeric_vector.h"
#include "quadrature.h"
#include "system.h"
#include "uniform_refinement_estimator.h"

#ifdef ENABLE_AMR

//-----------------------------------------------------------------
// ErrorEstimator implementations
void UniformRefinementEstimator::estimate_error (const System& _system,
					         ErrorVector& error_per_cell,
					         bool estimate_parent_error)
{
  START_LOG("estimate_error()", "UniformRefinementEstimator");
  this->_estimate_error
    (NULL, &_system, &error_per_cell, NULL, NULL, estimate_parent_error);
  STOP_LOG("estimate_error()", "UniformRefinementEstimator");
}

void UniformRefinementEstimator::estimate_errors (const EquationSystems& _es,
					         ErrorVector& error_per_cell,
                                                 std::map<const System*, std::vector<float> >& component_scales,
					         bool estimate_parent_error)
{
  START_LOG("estimate_errors()", "UniformRefinementEstimator");
  this->_estimate_error
    (&_es, NULL, &error_per_cell, NULL, &component_scales, estimate_parent_error);
  STOP_LOG("estimate_errors()", "UniformRefinementEstimator");
}

void UniformRefinementEstimator::estimate_errors (const EquationSystems& _es,
                                                  ErrorMap& errors_per_cell,
                                                  bool estimate_parent_error)
{
  START_LOG("estimate_errors()", "UniformRefinementEstimator");
  this->_estimate_error
    (&_es, NULL, NULL, &errors_per_cell, NULL, estimate_parent_error);
  STOP_LOG("estimate_errors()", "UniformRefinementEstimator");
}

void UniformRefinementEstimator::_estimate_error (const EquationSystems* _es,
                                                 const System* _system,
					         ErrorVector* error_per_cell,
                                                 ErrorMap* errors_per_cell,
                                                 std::map<const System*, std::vector<float> > *_component_scales,
					         bool)
{
  // Get a vector of the Systems we're going to work on,
  // and set up a component_scales map if necessary
  std::vector<System *> system_list;
  std::map<const System*, std::vector<float> > *component_scales;

  if (_es)
    {
      libmesh_assert(!_system);
      libmesh_assert(_es->n_systems());
      _system = &(_es->get_system(0));
      libmesh_assert(&(_system->get_equation_systems()) == _es);

      libmesh_assert(_es->n_systems());
      for (unsigned int i=0; i != _es->n_systems(); ++i)
      // We have to break the rules here, because we can't refine a const System
        system_list.push_back(const_cast<System *>(&(_es->get_system(i))));

      // If we're computing one vector, we need to know how to scale
      // each variable's contributions to it.
      if (_component_scales)
        {
          libmesh_assert(!errors_per_cell);
          component_scales = _component_scales;
        }
      else
      // If we're computing many vectors, we just need to know which
      // variables to skip
        {
          libmesh_assert (errors_per_cell);

          component_scales = new std::map<const System*, std::vector<float> >;
          for (unsigned int i=0; i!= _es->n_systems(); ++i)
            {
              const System &sys = _es->get_system(i);
              unsigned int n_vars = sys.n_vars();

              std::vector<float> cs(sys.n_vars(), 0.0);
              for (unsigned int v = 0; v != n_vars; ++v)
                {
                  if (errors_per_cell->find(std::make_pair(&sys, v)) ==
                      errors_per_cell->end())
                    continue;

                  cs[v] = 1.0;
                }
              (*component_scales)[&sys] = cs;
            }
        }
    }
  else
    {
      libmesh_assert(_system);
      // We have to break the rules here, because we can't refine a const System
      system_list.push_back(const_cast<System *>(_system));

      libmesh_assert(!_component_scales);
      component_scales = new std::map<const System*, std::vector<float> >;
      (*component_scales)[_system] = component_scale;
    }

  // An EquationSystems reference will be convenient.
  // We have to break the rules here, because we can't refine a const System
  EquationSystems& es =
    const_cast<EquationSystems &>(_system->get_equation_systems());

  // The current mesh
  MeshBase& mesh = es.get_mesh();

  // The dimensionality of the mesh
  const unsigned int dim = mesh.mesh_dimension();
  
  // Resize the error_per_cell vectors to be
  // the number of elements, initialize them to 0.
  if (error_per_cell)
    {
      error_per_cell->clear();
      error_per_cell->resize (mesh.max_elem_id(), 0.);
    }
  else
    {
      libmesh_assert(errors_per_cell);
      for (ErrorMap::iterator i = errors_per_cell->begin();
           i != errors_per_cell->end(); ++i)
        {
          ErrorVector *e = i->second;
          e->clear();
          e->resize(mesh.max_elem_id(), 0.);
        }
    }

  // component_mask has long since been deprecated
  if (!component_mask.empty())
    deprecated();

  // We'll want to back up all coarse grid vectors
  std::vector<std::map<std::string, NumericVector<Number> *> >
    coarse_vectors(system_list.size());
  std::vector<NumericVector<Number> *>
    coarse_solutions(system_list.size());
  std::vector<NumericVector<Number> *>
    coarse_local_solutions(system_list.size());
  // And make copies of projected solutions
  std::vector<NumericVector<Number> *>
    projected_solutions(system_list.size());

  // And we'll need to temporarily change solution projection settings
  std::vector<bool> old_projection_settings(system_list.size());

  for (unsigned int i=0; i != system_list.size(); ++i)
    {
      System &system = *system_list[i];


      // Check for valid component_scales
      if (!(*component_scales)[&system].empty())
        {
          if ((*component_scales)[&system].size() != system.n_vars())
	    {
	      std::cerr << "ERROR: component_scale is the wrong size:"
		        << std::endl
		        << " component_scales[" << i << "].scale()=" 
                        << (*component_scales)[&system].size()
		        << std::endl
		        << " n_vars=" << system.n_vars()
		        << std::endl;
	      libmesh_error();
	    }
        }
      else
        {
          // No specified scaling.  Scale all variables by one.
          (*component_scales)[&system].resize (system.n_vars());
          std::fill ((*component_scales)[&system].begin(),
                     (*component_scales)[&system].end(), 1.0);
        }
  
      // Back up the solution vector
      coarse_solutions[i] = system.solution->clone().release();
      coarse_local_solutions[i] =
        system.current_local_solution->clone().release();

      // Back up all other coarse grid vectors
      for (System::vectors_iterator vec = system.vectors_begin(); vec !=
           system.vectors_end(); ++vec)
        {
          // The (string) name of this vector
          const std::string& var_name = vec->first;

          coarse_vectors[i][var_name] = vec->second->clone().release();
        }

      // Make sure the solution is projected when we refine the mesh
      old_projection_settings[i] = system.project_solution_on_reinit();
      system.project_solution_on_reinit() = true;
    }

  // Find the number of coarse mesh elements, to make it possible
  // to find correct coarse elem ids later
  const unsigned int max_coarse_elem_id = mesh.max_elem_id();
#ifndef NDEBUG
  // n_coarse_elem is only used in an assertion later so
  // avoid declaring it unless asserts are active.
  const unsigned int n_coarse_elem = mesh.n_elem();
#endif
  
  // Uniformly refine the mesh
  MeshRefinement mesh_refinement(mesh);

  libmesh_assert (number_h_refinements > 0 || number_p_refinements > 0);

  // FIXME: this may break if there is more than one System
  // on this mesh but estimate_error was still called instead of
  // estimate_errors
  for (unsigned int i = 0; i != number_h_refinements; ++i)
    {
      mesh_refinement.uniformly_refine(1);
      es.reinit();
    }
      
  for (unsigned int i = 0; i != number_p_refinements; ++i)
    {
      mesh_refinement.uniformly_p_refine(1);
      es.reinit();
    }

  for (unsigned int i=0; i != system_list.size(); ++i)
    {
      System &system = *system_list[i];
      
      // Copy the projected coarse grid solutions, which will be
      // overwritten by solve()
//      projected_solutions[i] = system.solution->clone().release();