system.C

Go to the documentation of this file.

// $Id: system.C 3546 2009-11-06 20:13:28Z roystgnr $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
  
// 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 "dof_map.h"
#include "equation_systems.h"
#include "libmesh_logging.h"
#include "mesh_base.h"
#include "numeric_vector.h"
#include "parameter_vector.h"
#include "qoi_set.h"
#include "string_to_enum.h"
#include "system.h"
#include "utility.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()),
  qoi                      (0),
  _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 LIBMESH_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 ()
{
  _variables.clear();
  
  _variable_numbers.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();
  
  //If no variables have been added to this system
  //don't do anything
  if(!this->n_vars())
    return;
 
  // Then call the user-provided intialization function
  this->user_initialization();
}



void System::init_data ()
{
  MeshBase &mesh = this->get_mesh();

  // Distribute the degrees of freedom on the mesh
  _dof_map->distribute_dofs (mesh);

#ifdef LIBMESH_ENABLE_AMR

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

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

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

  // And clean up the send_list before we first use it
  _dof_map->prepare_send_list();

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

  // Resize the current_local_solution for the current mesh
#ifdef LIBMESH_ENABLE_GHOSTED
  current_local_solution->init (this->n_dofs(), this->n_local_dofs(),
                                _dof_map->get_send_list(), false,
                                GHOSTED);
#else
  current_local_solution->init (this->n_dofs(), false, SERIAL);
#endif

  // 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(), false, PARALLEL);
}



void System::restrict_vectors ()
{
#ifdef LIBMESH_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(), false, PARALLEL);
    }

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

  // Restrict the solution on the coarsened cells
  if (_solution_projection)
    this->project_vector (*solution);

#ifdef LIBMESH_ENABLE_GHOSTED
  current_local_solution->init(this->n_dofs(),
                               this->n_local_dofs(), send_list, 
                               false, GHOSTED);
#else
  current_local_solution->init(this->n_dofs());
#endif

  if (_solution_projection)
    solution->localize (*current_local_solution, send_list); 

#endif // LIBMESH_ENABLE_AMR
}



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



void System::reinit ()
{
#ifdef LIBMESH_ENABLE_AMR
  //If no variables have been added to this system
  //don't do anything
  if(!this->n_vars())
    return;

  // Constraints get handled in EquationSystems::reinit now
//  _dof_map->create_dof_constraints(this->get_mesh());

  // Update the solution based on the projected
  // current_local_solution.
  solution->init (this->n_dofs(), this->n_local_dofs(), true, PARALLEL);
        
  libmesh_assert (solution->size() == current_local_solution->size());
  // Not true with ghosted vectors
  // 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->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 ();

  // If this system is empty... don't do anything!
  if(!this->n_vars())
    return;

  // 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->size()       == solution->size());
  // Not true with ghosted vectors
  // libmesh_assert (current_local_solution->local_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");
}



void System::assemble_qoi (const QoISet& qoi_indices)
{
  // Log how long the user's matrix assembly code takes
  START_LOG("assemble_qoi()", "System");
  
  // Call the user-specified quantity of interest function
  this->user_QOI(qoi_indices);

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



void System::assemble_qoi_derivative (const QoISet& qoi_indices)
{
  // Log how long the user's matrix assembly code takes
  START_LOG("assemble_qoi_derivative()", "System");
  
  // Call the user-specified quantity of interest function
  this->user_QOI_derivative(qoi_indices);

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



void System::qoi_parameter_sensitivity
  (const QoISet& qoi_indices,
   const ParameterVector& parameters,
   SensitivityData& sensitivities)
{
  // Forward sensitivities are more efficient for Nq > Np
  if (qoi_indices.size(*this) > parameters.size())
    forward_qoi_parameter_sensitivity(qoi_indices, parameters, sensitivities);
  // Adjoint sensitivities are more efficient for Np > Nq,
  // and an adjoint may be more reusable than a forward 
  // solution sensitivity in the Np == Nq case.
  else
    adjoint_qoi_parameter_sensitivity(qoi_indices, parameters, sensitivities);
}



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;
        }
    }
  else
    overall_result = false;

  if (verbose)
    {
      std::cout << "   finished comparisons, ";
      if (overall_result)
        std::cout << "found no differences." << std::endl << std::endl;
      else 
        std::cout << "found differences." << std::endl << std::endl;
    }
          
  return overall_result;
}



void System::update_global_solution (std::vector<Number>& global_soln) const
{
  global_soln.resize (solution->size());

  solution->localize (global_soln);
}



void System::update_global_solution (std::vector<Number>& global_soln,
                                     const unsigned int   dest_proc) const
{
  global_soln.resize        (solution->size());

  solution->localize_to_one (global_soln, dest_proc);
}



NumericVector<Number> & System::add_vector (const std::string& vec_name,
                                            const bool projections)
{
  // Return the vector if it is already there.
  if (this->have_vector(vec_name))
    return *(_vectors[vec_name]);

  // Otherwise build the vector
  NumericVector<Number>* buf = NumericVector<Number>::build().release();
  _vectors.insert (std::make_pair (vec_name, buf));
  _vector_projections.insert (std::make_pair (vec_name, projections));

  // Initialize it if necessary
  if (!_can_add_vectors)
    buf->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);

  return *buf;
}



const NumericVector<Number> * System::request_vector (const std::string& vec_name) const
{
  const_vectors_iterator pos = _vectors.find(vec_name);
  
  if (pos == _vectors.end())
    return NULL;
  
  return pos->second;
}



NumericVector<Number> * System::request_vector (const std::string& vec_name)
{
  vectors_iterator pos = _vectors.find(vec_name);
  
  if (pos == _vectors.end())
    return NULL;
  
  return pos->second;
}



const NumericVector<Number> * System::request_vector (const unsigned int vec_num) const
{
  const_vectors_iterator v = vectors_begin();
  const_vectors_iterator v_end = vectors_end();
  unsigned int num = 0;
  while((num<vec_num) && (v!=v_end))
    {
      num++;
      ++v;
    }
  if (v==v_end)
    return NULL;
  return v->second;
}



NumericVector<Number> * System::request_vector (const unsigned int vec_num)
{
  vectors_iterator v = vectors_begin();
  vectors_iterator v_end = vectors_end();
  unsigned int num = 0;
  while((num<vec_num) && (v!=v_end))
    {
      num++;
      ++v;
    }
  if (v==v_end)
    return NULL;
  return v->second;
}



const NumericVector<Number> & System::get_vector (const std::string& vec_name) const
{
  // Make sure the vector exists
  const_vectors_iterator pos = _vectors.find(vec_name);
  
  if (pos == _vectors.end())
    {
      std::cerr << "ERROR: vector "
                << vec_name
                << " does not exist in this system!"
                << std::endl;      
      libmesh_error();
    }
  
  return *(pos->second);
}



NumericVector<Number> & System::get_vector (const std::string& vec_name)
{
  // Make sure the vector exists
  vectors_iterator pos = _vectors.find(vec_name);
  
  if (pos == _vectors.end())
    {
      std::cerr << "ERROR: vector "
                << vec_name
                << " does not exist in this system!"
                << std::endl;      
      libmesh_error();
    }
  
  return *(pos->second);
}



const NumericVector<Number> & System::get_vector (const unsigned int vec_num) const
{
  const_vectors_iterator v = vectors_begin();
  const_vectors_iterator v_end = vectors_end();
  unsigned int num = 0;
  while((num<vec_num) && (v!=v_end))
    {
      num++;
      ++v;
    }
  libmesh_assert(v!=v_end);
  return *(v->second);
}



NumericVector<Number> & System::get_vector (const unsigned int vec_num)
{
  vectors_iterator v = vectors_begin();
  vectors_iterator v_end = vectors_end();
  unsigned int num = 0;
  while((num<vec_num) && (v!=v_end))
    {
      num++;
      ++v;
    }
  libmesh_assert(v!=v_end);
  return *(v->second);
}



const std::string& System::vector_name (const unsigned int vec_num)
{
  const_vectors_iterator v = vectors_begin();
  const_vectors_iterator v_end = vectors_end();
  unsigned int num = 0;
  while((num<vec_num) && (v!=v_end))
    {
      num++;
      ++v;
    }
  libmesh_assert(v!=v_end);
  return v->first;
}



NumericVector<Number> & System::add_sensitivity_solution (unsigned int i)
{
  OStringStream sensitivity_name;
  sensitivity_name << "sensitivity_solution" << i;

  return this->add_vector(sensitivity_name.str());
}


NumericVector<Number> & System::get_sensitivity_solution (unsigned int i)
{
  OStringStream sensitivity_name;
  sensitivity_name << "sensitivity_solution" << i;

  return this->get_vector(sensitivity_name.str());
}



const NumericVector<Number> & System::get_sensitivity_solution (unsigned int i) const
{
  OStringStream sensitivity_name;
  sensitivity_name << "sensitivity_solution" << i;

  return this->get_vector(sensitivity_name.str());
}



NumericVector<Number> & System::add_weighted_sensitivity_solution ()
{
  return this->add_vector("weighted_sensitivity_solution");
}



NumericVector<Number> & System::get_weighted_sensitivity_solution ()
{
  return this->get_vector("weighted_sensitivity_solution");
}



const NumericVector<Number> & System::get_weighted_sensitivity_solution () const
{
  return this->get_vector("weighted_sensitivity_solution");
}



NumericVector<Number> & System::add_adjoint_solution (unsigned int i)
{
  OStringStream adjoint_name;
  adjoint_name << "adjoint_solution" << i;

  return this->add_vector(adjoint_name.str());
}



NumericVector<Number> & System::get_adjoint_solution (unsigned int i)
{
  OStringStream adjoint_name;
  adjoint_name << "adjoint_solution" << i;

  return this->get_vector(adjoint_name.str());
}



const NumericVector<Number> & System::get_adjoint_solution (unsigned int i) const
{
  OStringStream adjoint_name;
  adjoint_name << "adjoint_solution" << i;

  return this->get_vector(adjoint_name.str());
}



NumericVector<Number> & System::add_weighted_sensitivity_adjoint_solution (unsigned int i)
{
  OStringStream adjoint_name;
  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;

  return this->get_vector(adjoint_name.str());
}



NumericVector<Number> & System::get_weighted_sensitivity_adjoint_solution (unsigned int i)
{
  OStringStream adjoint_name;
  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;

  return this->get_vector(adjoint_name.str());
}



const NumericVector<Number> & System::get_weighted_sensitivity_adjoint_solution (unsigned int i) const
{
  OStringStream adjoint_name;
  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;

  return this->get_vector(adjoint_name.str());
}



NumericVector<Number> & System::add_adjoint_rhs (unsigned int i)
{
  OStringStream adjoint_rhs_name;
  adjoint_rhs_name << "adjoint_rhs" << i;

  return this->add_vector(adjoint_rhs_name.str());
}



NumericVector<Number> & System::get_adjoint_rhs (unsigned int i)
{
  OStringStream adjoint_rhs_name;
  adjoint_rhs_name << "adjoint_rhs" << i;

  return this->get_vector(adjoint_rhs_name.str());
}



const NumericVector<Number> & System::get_adjoint_rhs (unsigned int i) const
{
  OStringStream adjoint_rhs_name;
  adjoint_rhs_name << "adjoint_rhs" << i;

  return this->get_vector(adjoint_rhs_name.str());
}



NumericVector<Number> & System::add_sensitivity_rhs (unsigned int i)
{
  OStringStream sensitivity_rhs_name;
  sensitivity_rhs_name << "sensitivity_rhs" << i;

  return this->add_vector(sensitivity_rhs_name.str());
}



NumericVector<Number> & System::get_sensitivity_rhs (unsigned int i)
{
  OStringStream sensitivity_rhs_name;
  sensitivity_rhs_name << "sensitivity_rhs" << i;

  return this->get_vector(sensitivity_rhs_name.str());
}



const NumericVector<Number> & System::get_sensitivity_rhs (unsigned int i) const
{
  OStringStream sensitivity_rhs_name;
  sensitivity_rhs_name << "sensitivity_rhs" << i;

  return this->get_vector(sensitivity_rhs_name.str());
}



unsigned int System::add_variable (const std::string& var,
                                   const FEType& type,
                                   const std::set<subdomain_id_type> * const active_subdomains)
{  
  // Make sure the variable isn't there already
  // or if it is, that it's the type we want
  for (unsigned int v=0; v<this->n_vars(); v++)
    if (this->variable_name(v) == var)
      {
        if (this->variable_type(v) == type)
          return _variables[v].number();

        std::cerr << "ERROR: incompatible variable "
                  << var
                  << " has already been added for this system!"
                  << std::endl;
        libmesh_error();
      }

  const unsigned int curr_n_vars = this->n_vars();                                              

  // Add the variable to the list
  _variables.push_back((active_subdomains == NULL) ?
                       Variable(var, curr_n_vars, type) :
                       Variable(var, curr_n_vars, type, *active_subdomains));

  libmesh_assert ((curr_n_vars+1) == this->n_vars());

  _variable_numbers[var] = curr_n_vars;

  // Add the variable to the _dof_map
  _dof_map->add_variable (_variables.back());

  // Return the number of the new variable
  return curr_n_vars;
}



unsigned int System::add_variable (const std::string& var,
                                   const Order order,
                                   const FEFamily family,
                                   const std::set<subdomain_id_type> * const active_subdomains)
{
  return this->add_variable(var, 
                            FEType(order, family), 
                            active_subdomains);
}



bool System::has_variable (const std::string& var) const
{
  return _variable_numbers.count(var);
}



unsigned short int System::variable_number (const std::string& var) const
{
  // Make sure the variable exists
  std::map<std::string, unsigned short int>::const_iterator
    pos = _variable_numbers.find(var);
  
  if (pos == _variable_numbers.end())
    {
      std::cerr << "ERROR: variable "
                << var
                << " does not exist in this system!"
                << std::endl;      
      libmesh_error();
    }
  libmesh_assert (_variables[pos->second].name() == var);

  return pos->second;
}


void System::local_dof_indices(const unsigned int var, std::set<unsigned int> & var_indices) const
{
  // Make sure the set is clear
  var_indices.clear();

  std::vector<unsigned int> dof_indices;
  
  // Begin the loop over the elements
  MeshBase::const_element_iterator       el     =
    this->get_mesh().active_local_elements_begin();
  const MeshBase::const_element_iterator end_el =
    this->get_mesh().active_local_elements_end();

  const unsigned int 
    first_local = this->get_dof_map().first_dof(),
    end_local   = this->get_dof_map().end_dof();

  for ( ; el != end_el; ++el)
    {
      const Elem* elem = *el;      
      this->get_dof_map().dof_indices (elem, dof_indices, var);

      for(unsigned int i=0; i<dof_indices.size(); i++)
        {
          unsigned int dof = dof_indices[i];

          //If the dof is owned by the local processor
          if(first_local <= dof && dof < end_local)
            var_indices.insert(dof_indices[i]);
        }
    }
}

void System::zero_variable (NumericVector<Number>& v, unsigned int var_num) const
{
  /* Make sure the call makes sense.  */
  libmesh_assert(var_num<this->n_vars());

  /* Get a reference to the mesh.  */
  const MeshBase& mesh = this->get_mesh();

  /* Check which system we are.  */
  const unsigned int sys_num = this->number();

  /* Loop over nodes.  */
  {
    MeshBase::const_node_iterator it = mesh.local_nodes_begin();
    const MeshBase::const_node_iterator end_it = mesh.local_nodes_end(); 
    for ( ; it != end_it; ++it)
      {    
        const Node* node = *it;
        unsigned int n_comp = node->n_comp(sys_num,var_num);
        for(unsigned int i=0; i<n_comp; i++)
          {
            const unsigned int index = node->dof_number(sys_num,var_num,i);
            v.set(index,0.0);
          }
      }
  }

  /* Loop over elements.  */
  {
    MeshBase::const_element_iterator it = mesh.active_local_elements_begin();
    const MeshBase::const_element_iterator end_it = mesh.active_local_elements_end(); 
    for ( ; it != end_it; ++it)
      {    
        const Elem* elem = *it;
        unsigned int n_comp = elem->n_comp(sys_num,var_num);
        for(unsigned int i=0; i<n_comp; i++)
          {
            const unsigned int index = elem->dof_number(sys_num,var_num,i);
            v.set(index,0.0);
          }
      }
  }
}

Real System::discrete_var_norm(NumericVector<Number>& v,
                               unsigned int var,
                               FEMNormType norm_type) const
{
  std::set<unsigned int> var_indices;
  local_dof_indices(var, var_indices);

  if(norm_type == DISCRETE_L1)
    return v.subset_l1_norm(var_indices);
  if(norm_type == DISCRETE_L2)
    return v.subset_l2_norm(var_indices);
  if(norm_type == DISCRETE_L_INF)
    return v.subset_linfty_norm(var_indices);
  else
    libmesh_error();
}

Real System::calculate_norm(NumericVector<Number>& v,
                            unsigned int var,
                            FEMNormType norm_type) const
{
  //short circuit to save time
  if(norm_type == DISCRETE_L1 ||
     norm_type == DISCRETE_L2 ||
     norm_type == DISCRETE_L_INF)
    return discrete_var_norm(v,var,norm_type);

  // Not a discrete norm
  std::vector<FEMNormType> norms(this->n_vars(), L2);
  std::vector<Real> weights(this->n_vars(), 0.0);
  norms[var] = norm_type;
  weights[var] = 1.0;
  Real val = this->calculate_norm(v, SystemNorm(norms, weights));
  return val;
}



Real System::calculate_norm(NumericVector<Number>& v,
                            const SystemNorm &norm) const
{
  // This function must be run on all processors at once
  parallel_only();

  START_LOG ("calculate_norm()", "System");

  // Zero the norm before summation
  Real v_norm = 0.;

  if (norm.is_discrete())
    {
      STOP_LOG ("calculate_norm()", "System");
      //Check to see if all weights are 1.0
      unsigned int check_var = 0;
      for (; check_var != this->n_vars(); ++check_var)
        if(norm.weight(check_var) != 1.0)
          break;

      //All weights were 1.0 so just do the full vector discrete norm
      if(check_var == this->n_vars())
        {
          FEMNormType norm_type = norm.type(0);
          
          if(norm_type == DISCRETE_L1)
            return v.l1_norm();
          if(norm_type == DISCRETE_L2)
            return v.l2_norm();
          if(norm_type == DISCRETE_L_INF)
            return v.linfty_norm();
          else
            libmesh_error();
        }

      for (unsigned int var=0; var != this->n_vars(); ++var)
        {
          // Skip any variables we don't need to integrate
          if (norm.weight(var) == 0.0)
            continue;

          v_norm += norm.weight(var) * discrete_var_norm(v, var, norm.type(var));
        }

      return v_norm;
    }

  // Localize the potentially parallel vector
  AutoPtr<NumericVector<Number> > local_v = NumericVector<Number>::build();
  local_v->init(v.size(), true, SERIAL);
  v.localize (*local_v, _dof_map->get_send_list()); 

  unsigned int dim = this->get_mesh().mesh_dimension();

  // Loop over all variables
  for (unsigned int var=0; var != this->n_vars(); ++var)
    {
      // Skip any variables we don't need to integrate
      if (norm.weight(var) == 0.0)
        continue;

      const FEType& fe_type = this->get_dof_map().variable_type(var);
      AutoPtr<QBase> qrule =
        fe_type.default_quadrature_rule (dim);
      AutoPtr<FEBase> fe
        (FEBase::build(dim, fe_type));
      fe->attach_quadrature_rule (qrule.get());

      const std::vector<Real>&               JxW = fe->get_JxW();
      const std::vector<std::vector<Real> >* phi = NULL;
      if (norm.type(var) == H1 ||
          norm.type(var) == H2 ||
          norm.type(var) == L2)
        phi = &(fe->get_phi());

      const std::vector<std::vector<RealGradient> >* dphi = NULL;
      if (norm.type(var) == H1 ||
          norm.type(var) == H2 ||
          norm.type(var) == H1_SEMINORM)
        dphi = &(fe->get_dphi());
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
      const std::vector<std::vector<RealTensor> >*   d2phi = NULL;
      if (norm.type(var) == H2 ||
          norm.type(var) == H2_SEMINORM)
        d2phi = &(fe->get_d2phi());
#endif

      std::vector<unsigned int> dof_indices;

      // Begin the loop over the elements
      MeshBase::const_element_iterator       el     =
        this->get_mesh().active_local_elements_begin();
      const MeshBase::const_element_iterator end_el =
        this->get_mesh().active_local_elements_end();

      for ( ; el != end_el; ++el)
        {
          const Elem* elem = *el;

          fe->reinit (elem);

          this->get_dof_map().dof_indices (elem, dof_indices, var);

          const unsigned int n_qp = qrule->n_points();

          const unsigned int n_sf = dof_indices.size();

          // Begin the loop over the Quadrature points.
          for (unsigned int qp=0; qp<n_qp; qp++)
            {
              if (norm.type(var) == H1 ||
                  norm.type(var) == H2 ||
                  norm.type(var) == L2)
                {
                  Number u_h = 0.;
                  for (unsigned int i=0; i != n_sf; ++i)
                    u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
                  v_norm += norm.weight_sq(var) *
                            JxW[qp] * libmesh_norm(u_h);
                }

              if (norm.type(var) == H1 ||
                  norm.type(var) == H2 ||
                  norm.type(var) == H1_SEMINORM)
                {
                  Gradient grad_u_h;
                  for (unsigned int i=0; i != n_sf; ++i)
                    grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
                  v_norm += norm.weight_sq(var) *
                            JxW[qp] * grad_u_h.size_sq();
                }

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
              if (norm.type(var) == H2 ||
                  norm.type(var) == H2_SEMINORM)
                {
                  Tensor hess_u_h;
                  for (unsigned int i=0; i != n_sf; ++i)
                    hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
                  v_norm += norm.weight_sq(var) *
                            JxW[qp] * hess_u_h.size_sq();
                }
#endif
            }
        }
    }

  Parallel::sum(v_norm);

  STOP_LOG ("calculate_norm()", "System");

  return std::sqrt(v_norm);
}



std::string System::get_info() const
{
  std::ostringstream out;

  
  const std::string& sys_name = this->name();
      
  out << "   System \"" << sys_name << "\"\n"
      << "    Type \""  << this->system_type() << "\"\n"
      << "    Variables=";
  
  for (unsigned int vn=0; vn<this->n_vars(); vn++)
      out << "\"" << this->variable_name(vn) << "\" ";
     
  out << '\n';

  out << "    Finite Element Types=";
#ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
  for (unsigned int vn=0; vn<this->n_vars(); vn++)
    out << "\""
        << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_type(vn).family)
        << "\" ";  
#else
  for (unsigned int vn=0; vn<this->n_vars(); vn++)
    {
      out << "\""
          << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_type(vn).family)
          << "\", \""
          << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_type(vn).radial_family)
          << "\" ";
    }

  out << '\n' << "    Infinite Element Mapping=";
  for (unsigned int vn=0; vn<this->n_vars(); vn++)
    out << "\""
        << Utility::enum_to_string<InfMapType>(this->get_dof_map().variable_type(vn).inf_map)
        << "\" ";
#endif      

  out << '\n';
      
  out << "    Approximation Orders=";
  for (unsigned int vn=0; vn<this->n_vars(); vn++)
    {
#ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
      out << "\""
          << Utility::enum_to_string<Order>(this->get_dof_map().variable_type(vn).order)
          << "\" ";
#else
      out << "\""
          << Utility::enum_to_string<Order>(this->get_dof_map().variable_type(vn).order)
          << "\", \""
          << Utility::enum_to_string<Order>(this->get_dof_map().variable_type(vn).radial_order)
          << "\" ";
#endif
    }

  out << '\n';
      
  out << "    n_dofs()="             << this->n_dofs()             << '\n';
  out << "    n_local_dofs()="       << this->n_local_dofs()       << '\n';
#ifdef LIBMESH_ENABLE_AMR
  out << "    n_constrained_dofs()=" << this->n_constrained_dofs() << '\n';
#endif

  out << "    " << "n_vectors()="  << this->n_vectors()  << '\n';
//   out << "    " << "n_additional_matrices()=" << this->n_additional_matrices() << '\n';
  
  return out.str();
}



void System::attach_init_function (void fptr(EquationSystems& es,
                                             const std::string& name))
{
  libmesh_assert (fptr != NULL);
  
  _init_system = fptr;
}



void System::attach_assemble_function (void fptr(EquationSystems& es,
                                                 const std::string& name))
{
  libmesh_assert (fptr != NULL);
  
  _assemble_system = fptr;  
}



void System::attach_constraint_function(void fptr(EquationSystems& es,
                                                  const std::string& name))
{
  libmesh_assert (fptr != NULL);
  
  _constrain_system = fptr;  
}



void System::attach_QOI_function(void fptr(EquationSystems&,
                                           const std::string&,
                                           const QoISet&))
{
  libmesh_assert (fptr != NULL);
  
  _qoi_evaluate = fptr;  
}



void System::attach_QOI_derivative(void fptr(EquationSystems&,
                                             const std::string&,
                                             const QoISet&))
{
  libmesh_assert (fptr != NULL);
  
  _qoi_evaluate_derivative = fptr;  
}



void System::user_initialization ()
{
  // Call the user-provided intialization function,
  // if it was provided
  if (_init_system != NULL)
    this->_init_system (_equation_systems, this->name());
}


void System::user_assembly ()
{
  // Call the user-provided assembly function,
  // if it was provided
  if (_assemble_system != NULL)
    this->_assemble_system (_equation_systems, this->name());
}



void System::user_constrain ()
{
  // Call the user-provided constraint function, 
  // if it was provided
  if(_constrain_system!= NULL)
    this->_constrain_system(_equation_systems, this->name());
}



void System::user_QOI (const QoISet& qoi_indices)
{
  // Call the user-provided quantity of interest function, 
  // if it was provided
  if(_qoi_evaluate != NULL)
    this->_qoi_evaluate(_equation_systems, this->name(), qoi_indices);
}



void System::user_QOI_derivative (const QoISet& qoi_indices)
{
  // Call the user-provided quantity of interest derivative, 
  // if it was provided
  if(_qoi_evaluate_derivative != NULL)
    this->_qoi_evaluate_derivative(_equation_systems, this->name(), qoi_indices);
}