system.c

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

// $Id: system.C 2983 2008-08-15 20:36:36Z 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 <sstream>   // for std::ostringstream


// Local includes
#include "system.h"
#include "equation_systems.h"
#include "libmesh_logging.h"
#include "utility.h"
#include "string_to_enum.h"
#include "dof_map.h"
#include "numeric_vector.h"
#include "mesh_base.h"

// includes for calculate_norm
#include "fe_base.h"
#include "parallel.h"
#include "quadrature.h"
#include "tensor_value.h"
#include "vector_value.h"


// ------------------------------------------------------------
// System implementation
System::System (EquationSystems& es,
		const std::string& name,
		const unsigned int number) :

  assemble_before_solve    (true),
  solution                 (NumericVector<Number>::build()),
  current_local_solution   (NumericVector<Number>::build()),
  _init_system             (NULL),
  _assemble_system         (NULL),
  _constrain_system        (NULL),
  _dof_map                 (new DofMap(number)),
  _equation_systems        (es),
  _mesh                    (es.get_mesh()),
  _sys_name                (name),
  _sys_number              (number),
  _active                  (true),
  _solution_projection     (true),
  _can_add_vectors         (true),
  _additional_data_written (false)
{
}



System::~System ()
{
  // Null-out the function pointers.  Since this
  // class is getting destructed it is pointless,
  // but a good habit.
  _init_system = _assemble_system = NULL;

  // Clear data
  this->clear ();

  libmesh_assert (!libMesh::closed());
}




unsigned int System::n_dofs() const
{
  return _dof_map->n_dofs();
}



unsigned int System::n_constrained_dofs() const
{
#ifdef ENABLE_AMR

  return _dof_map->n_constrained_dofs();

#else

  return 0;

#endif
}





unsigned int System::n_local_dofs() const
{
  return _dof_map->n_dofs_on_processor (libMesh::processor_id());
}






Number System::current_solution (const unsigned int global_dof_number) const
{
  // Check the sizes
  libmesh_assert (global_dof_number < _dof_map->n_dofs());
  libmesh_assert (global_dof_number < current_local_solution->size());
   
  return (*current_local_solution)(global_dof_number);
}





void System::clear ()
{
  _var_names.clear ();

  _var_type.clear ();
  
  _var_num.clear ();
  
  _dof_map->clear ();
  
  solution->clear ();

  current_local_solution->clear ();
  
  // clear any user-added vectors
  {
    for (vectors_iterator pos = _vectors.begin(); pos != _vectors.end(); ++pos)
      {
	pos->second->clear ();
	delete pos->second;
	pos->second = NULL;
      }
    
    _vectors.clear();
    _vector_projections.clear();
    _can_add_vectors = true;
  }

}



void System::init ()
{
  // First initialize any required data
  this->init_data();

  // Then call the user-provided intialization function
  this->user_initialization();
}



void System::init_data ()
{
  // Distribute the degrees of freedom on the mesh
  _dof_map->distribute_dofs (this->get_mesh());

#ifdef ENABLE_AMR

  // Recreate any hanging node constraints
  _dof_map->create_dof_constraints(this->get_mesh());

  // Apply any user-defined constraints
  this->user_constrain();

  // Expand any recursive constraints
  _dof_map->process_recursive_constraints();

#endif
  
  // Resize the solution conformal to the current mesh
  solution->init (this->n_dofs(), this->n_local_dofs());

  // Resize the current_local_solution for the current mesh
  current_local_solution->init (this->n_dofs());

  // from now on, no chance to add additional vectors
  _can_add_vectors = false;

  // initialize & zero other vectors, if necessary
  for (vectors_iterator pos = _vectors.begin(); pos != _vectors.end(); ++pos)
      pos->second->init (this->n_dofs(), this->n_local_dofs());
}



void System::restrict_vectors ()
{
#ifdef ENABLE_AMR
  // Restrict the _vectors on the coarsened cells
  for (vectors_iterator pos = _vectors.begin(); pos != _vectors.end(); ++pos)
    {
      NumericVector<Number>* v = pos->second;
      
      if (_vector_projections[pos->first])
	this->project_vector (*v);
      else
        v->init (this->n_dofs(), this->n_local_dofs());
    }

  // Restrict the solution on the coarsened cells
  if (_solution_projection)
    {
      this->project_vector (*solution);
      current_local_solution->clear();
      current_local_solution->init(this->n_dofs());
      const std::vector<unsigned int>& send_list = _dof_map->get_send_list ();
      solution->localize (*current_local_solution, send_list); 
    }
  else
    {
      current_local_solution->clear();
      current_local_solution->init(this->n_dofs());
    }
#endif
}



void System::prolong_vectors ()
{
#ifdef ENABLE_AMR
  // Currently project_vector handles both restriction and prolongation
  this->restrict_vectors();
#endif
}



void System::reinit ()
{
#ifdef ENABLE_AMR
  // Recreate any hanging node constraints
//  _dof_map->create_dof_constraints(this->get_mesh());

  // Apply any user-defined constraints
  this->user_constrain();

  // Expand any recursive constraints
  _dof_map->process_recursive_constraints();

  // Update the solution based on the projected
  // current_local_solution.
  solution->init (this->n_dofs(), this->n_local_dofs());
	
  libmesh_assert (solution->size() == current_local_solution->size());
  libmesh_assert (solution->size() == current_local_solution->local_size());

  const unsigned int first_local_dof = solution->first_local_index();
  const unsigned int local_size      = solution->local_size();
      
  for (unsigned int i=0; i<local_size; i++)
    solution->set(i+first_local_dof,
		  (*current_local_solution)(i+first_local_dof));
#endif
}



void System::update ()
{
  const std::vector<unsigned int>& send_list = _dof_map->get_send_list ();

  // Check sizes
  libmesh_assert (current_local_solution->local_size() == solution->size());
// More processors than elements => empty send_list
//  libmesh_assert (!send_list.empty());
  libmesh_assert (send_list.size() <= solution->size());

  // Create current_local_solution from solution.  This will
  // put a local copy of solution into current_local_solution.
  // Only the necessary values (specified by the send_list)
  // are copied to minimize communication
  solution->localize (*current_local_solution, send_list); 
}



void System::re_update ()
{
  //const std::vector<unsigned int>& send_list = _dof_map->get_send_list ();

  // Explicitly build a send_list
  std::vector<unsigned int> send_list(solution->size());
  Utility::iota (send_list.begin(), send_list.end(), 0);
  
  // Check sizes
  libmesh_assert (current_local_solution->local_size() == solution->size());
  libmesh_assert (current_local_solution->size()       == solution->size());
  libmesh_assert (!send_list.empty());
  libmesh_assert (send_list.size() <= solution->size());

  // Create current_local_solution from solution.  This will
  // put a local copy of solution into current_local_solution.
  solution->localize (*current_local_solution, send_list);
}



void System::assemble ()
{
  // Log how long the user's matrix assembly code takes
  START_LOG("assemble()", "System");
  
  // Call the user-specified assembly function
  this->user_assembly();

  // Stop logging the user code
  STOP_LOG("assemble()", "System");
}



bool System::compare (const System& other_system, 
		      const Real threshold,
		      const bool verbose) const
{
  // we do not care for matrices, but for vectors
  libmesh_assert (!_can_add_vectors);
  libmesh_assert (!other_system._can_add_vectors);

  if (verbose)
    {
      std::cout << "  Systems \"" << _sys_name << "\"" << std::endl;
      std::cout << "   comparing matrices not supported." << std::endl;
      std::cout << "   comparing names...";
    }

  // compare the name: 0 means identical
  const int name_result = _sys_name.compare(other_system.name());
  if (verbose)
    {
      if (name_result == 0)
	std::cout << " identical." << std::endl;
      else
	std::cout << "  names not identical." << std::endl;
      std::cout << "   comparing solution vector...";
    }


  // compare the solution: -1 means identical
  const int solu_result = solution->compare (*other_system.solution.get(),
					     threshold);

  if (verbose)
    {
      if (solu_result == -1)
	std::cout << " identical up to threshold." << std::endl;
      else
	std::cout << "  first difference occured at index = " 
		  << solu_result << "." << std::endl;
    }


  // safety check, whether we handle at least the same number
  // of vectors
  std::vector<int> ov_result;

  if (this->n_vectors() != other_system.n_vectors())
    {
      if (verbose)
        {
	  std::cout << "   Fatal difference. This system handles " 
		    << this->n_vectors() << " add'l vectors," << std::endl
		    << "   while the other system handles "
		    << other_system.n_vectors() 
		    << " add'l vectors." << std::endl
		    << "   Aborting comparison." << std::endl;
	}
      return false;
    }
  else if (this->n_vectors() == 0)
    {
      // there are no additional vectors...
      ov_result.clear ();
    }
  else
    {
      // compare other vectors
      for (const_vectors_iterator pos = _vectors.begin();
	   pos != _vectors.end(); ++pos)
        {
	  if (verbose)
	      std::cout << "   comparing vector \""
			<< pos->first << "\" ...";

	  // assume they have the same name
	  const NumericVector<Number>& other_system_vector = 
	      other_system.get_vector(pos->first);

	  ov_result.push_back(pos->second->compare (other_system_vector,
						    threshold));

	  if (verbose)
	    {
	      if (ov_result[ov_result.size()-1] == -1)
		std::cout << " identical up to threshold." << std::endl;
	      else
		std::cout << " first difference occured at" << std::endl
			  << "   index = " << ov_result[ov_result.size()-1] << "." << std::endl;
	    }

	}

    } // finished comparing additional vectors


  bool overall_result;
       
  // sum up the results
  if ((name_result==0) && (solu_result==-1))
    {
      if (ov_result.size()==0)
	overall_result = true;
      else
        {
	  bool ov_identical;
	  unsigned int n    = 0;
	  do
	    {
	      ov_identical = (ov_result[n]==-1);
	      n++;
	    }
	  while (ov_identical && n<ov_result.size());
	  overall_result = ov_identical;
	}