system.c

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

    }
  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]);
    }

  // We can only add new vectors before initializing...
  if (!_can_add_vectors)
    {
      std::cerr << "ERROR: Too late.  Cannot add vectors to the system after initialization"
		<< std::endl
		<< " any more.  You should have done this earlier."
		<< std::endl;
      libmesh_error();
    }

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

  return *buf;
}



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);
}



void System::add_variable (const std::string& var,
			   const FEType& type)
{  
  // Make sure the variable isn't there already
  // or if it is, that it's the type we want
  if (_var_num.count(var))
    {
      if (_var_type[var] == type)
        return;

      std::cerr << "ERROR: incompatible variable "
		<< var
		<< " has already been added for this system!"
		<< std::endl;
      libmesh_error();
    }
  
  // Add the variable to the list
  _var_names.push_back (var);
  _var_type[var]  = type;
  _var_num[var]   = (this->n_vars()-1);

  // Add the variable to the _dof_map
  _dof_map->add_variable (type);
}



void System::add_variable (const std::string& var,
			       const Order order,
			       const FEFamily family)
{
  FEType fe_type(order, family);
  
  this->add_variable(var, fe_type);
}



bool System::has_variable (const std::string& var) const
{
  if (_var_num.find(var) == _var_num.end())
    return false;
  return true;
}



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 = _var_num.find(var);
  
  if (pos == _var_num.end())
    {
      std::cerr << "ERROR: variable "
		<< var
		<< " does not exist in this system!"
		<< std::endl;      
      libmesh_error();
    }
  
  return pos->second;
}



Real System::calculate_norm(NumericVector<Number>& v,
                            unsigned int var,
                            FEMNormType norm_type) const
{
  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");

  if (norm.is_discrete())
    {
      STOP_LOG ("calculate_norm()", "System");
      if (norm.type(0) == DISCRETE_L1)
        return v.l1_norm();
      if (norm.type(0) == DISCRETE_L2)
        return v.l2_norm();
      if (norm.type(0) == DISCRETE_L_INF)
        return v.linfty_norm();
    }

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

  // Localize the potentially parallel vector
  AutoPtr<NumericVector<Number> > local_v = NumericVector<Number>::build();
  local_v->init(v.size(), v.size());
  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;
      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 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(var) * norm.weight(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(var) * norm.weight(var) *
                            JxW[qp] * grad_u_h.size_sq();
                }

#ifdef 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(var) * norm.weight(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 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 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 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::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());
}