equation_systems.c

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

//       const Elem* elem = *it;

//       // Fill the dof_indices vector for this variable
//       system.get_dof_map().dof_indices(elem,
// 				       dof_indices,
// 				       variable_num);

//       // Resize the element solution vector to fit the
//       // dof_indices for this element.
//       elem_soln.resize(dof_indices.size());

//       // Transfer the system solution to the element
//       // solution by mapping it through the dof_indices vector.
//       for (unsigned int i=0; i<dof_indices.size(); i++)
// 	elem_soln[i] = sys_soln[dof_indices[i]];

//       // Using the FE interface, compute the nodal_soln
//       // for the current elemnt type given the elem_soln
//       FEInterface::nodal_soln (dim,
// 			       fe_type,
// 			       elem,
// 			       elem_soln,
// 			       nodal_soln);

//       // Sanity check -- make sure that there are the same number
//       // of entries in the nodal_soln as there are nodes in the
//       // element!
//       libmesh_assert (nodal_soln.size() == elem->n_nodes());

//       // Copy the nodal solution over into the correct place in
//       // the global soln vector which will be returned to the user.
//       for (unsigned int n=0; n<elem->n_nodes(); n++)
// 	soln[elem->node(n)] = nodal_soln[n];
//     }
}




void EquationSystems::build_solution_vector (std::vector<Number>& soln) const
{
  // This function must be run on all processors at once
  parallel_only();

  libmesh_assert (this->n_systems());

  const unsigned int dim = _mesh.mesh_dimension();
  const unsigned int nn  = _mesh.n_nodes();
  const unsigned int nv  = this->n_vars();

  // We'd better have a contiguous node numbering
  libmesh_assert (nn == _mesh.max_node_id());

  // allocate storage to hold
  // (number_of_nodes)*(number_of_variables) entries.
  soln.resize(nn*nv);

  std::vector<Number>             sys_soln;

  // (Note that we use an unsigned short int here even though an
  // unsigned char would be more that sufficient.  The MPI 1.1
  // standard does not require that MPI_SUM, MPI_PROD etc... be
  // implemented for char data types. 12/23/2003 - BSK)  
  std::vector<unsigned short int> node_conn(nn);

  
  // Zero out the soln vector
  std::fill (soln.begin(),       soln.end(),       libMesh::zero);

  
  // Get the number of elements that share each node.  We will
  // compute the average value at each node.  This is particularly
  // useful for plotting discontinuous data.
  MeshBase::element_iterator       e_it  = _mesh.active_local_elements_begin();
  const MeshBase::element_iterator e_end = _mesh.active_local_elements_end(); 

  for ( ; e_it != e_end; ++e_it)
    for (unsigned int n=0; n<(*e_it)->n_nodes(); n++)
      node_conn[(*e_it)->node(n)]++;

  // Gather the distributed node_conn arrays in the case of
  // multiple processors
  Parallel::sum(node_conn);

  unsigned int var_num=0;

  // For each system in this EquationSystems object,
  // update the global solution and if we are on processor 0,
  // loop over the elements and build the nodal solution
  // from the element solution.  Then insert this nodal solution
  // into the vector passed to build_solution_vector.
  const_system_iterator       pos = _systems.begin();
  const const_system_iterator end = _systems.end();

  for (; pos != end; ++pos)
    {  
      const System& system  = *(pos->second);
      const unsigned int nv_sys = system.n_vars();
      
      system.update_global_solution (sys_soln);
      
      std::vector<Number>       elem_soln;   // The finite element solution
      std::vector<Number>       nodal_soln;  // The FE solution interpolated to the nodes
      std::vector<unsigned int> dof_indices; // The DOF indices for the finite element 
      
      for (unsigned int var=0; var<nv_sys; var++)
	{
	  const FEType& fe_type    = system.variable_type(var);
	  
	  MeshBase::element_iterator       it  = _mesh.active_local_elements_begin();
	  const MeshBase::element_iterator end = _mesh.active_local_elements_end(); 

	  for ( ; it != end; ++it)
	    {
	      const Elem* elem = *it;
	      system.get_dof_map().dof_indices (elem, dof_indices, var);
	      
	      elem_soln.resize(dof_indices.size());
	      
	      for (unsigned int i=0; i<dof_indices.size(); i++)
		elem_soln[i] = sys_soln[dof_indices[i]];
		  
	      FEInterface::nodal_soln (dim,
				       fe_type,
				       elem,
				       elem_soln,
				       nodal_soln);

#ifdef ENABLE_INFINITE_ELEMENTS
	      // infinite elements should be skipped...
	      if (!elem->infinite())
#endif
		{ 
		  libmesh_assert (nodal_soln.size() == elem->n_nodes());
		  
		  for (unsigned int n=0; n<elem->n_nodes(); n++)
		    {
		      libmesh_assert (node_conn[elem->node(n)] != 0);
		      soln[nv*(elem->node(n)) + (var + var_num)] +=
			nodal_soln[n]/static_cast<Real>(node_conn[elem->node(n)]);
		    }
		}
	    }	 
	}

      var_num += nv_sys;
    }

  // Now each processor has computed contriburions to the
  // soln vector.  Gather them all up.
  Parallel::sum(soln);
}




void EquationSystems::build_discontinuous_solution_vector (std::vector<Number>& soln) const
{
  libmesh_assert (this->n_systems());

  const unsigned int dim = _mesh.mesh_dimension();
  const unsigned int nv  = this->n_vars();
  unsigned int tw=0;

  // get the total weight
  {
    MeshBase::element_iterator       it  = _mesh.active_elements_begin();
    const MeshBase::element_iterator end = _mesh.active_elements_end(); 

    for ( ; it != end; ++it)
      tw += (*it)->n_nodes();
  }
  

  // Only if we are on processor zero, allocate the storage
  // to hold (number_of_nodes)*(number_of_variables) entries.
  if (_mesh.processor_id() == 0)
    soln.resize(tw*nv);

  std::vector<Number> sys_soln; 

  
  unsigned int var_num=0;

  // For each system in this EquationSystems object,
  // update the global solution and if we are on processor 0,
  // loop over the elements and build the nodal solution
  // from the element solution.  Then insert this nodal solution
  // into the vector passed to build_solution_vector.
  const_system_iterator       pos = _systems.begin();
  const const_system_iterator end = _systems.end();

  for (; pos != end; ++pos)
    {  
      const System& system  = *(pos->second);
      const unsigned int nv_sys = system.n_vars();
      
      system.update_global_solution (sys_soln, 0);
      
      if (_mesh.processor_id() == 0)
	{
	  std::vector<Number>       elem_soln;   // The finite element solution
	  std::vector<Number>       nodal_soln;  // The FE solution interpolated to the nodes
	  std::vector<unsigned int> dof_indices; // The DOF indices for the finite element 
	      
	  for (unsigned int var=0; var<nv_sys; var++)
	    {
	      const FEType& fe_type    = system.variable_type(var);

	      MeshBase::element_iterator       it  = _mesh.active_elements_begin();
	      const MeshBase::element_iterator end = _mesh.active_elements_end(); 

	      unsigned int nn=0;
	      
	      for ( ; it != end; ++it)
		{
		  const Elem* elem = *it;
		  system.get_dof_map().dof_indices (elem, dof_indices, var);
		  
		  elem_soln.resize(dof_indices.size());
		  
		  for (unsigned int i=0; i<dof_indices.size(); i++)
		    elem_soln[i] = sys_soln[dof_indices[i]];
		  
		  FEInterface::nodal_soln (dim,
					   fe_type,
					   elem,
					   elem_soln,
					   nodal_soln);

#ifdef ENABLE_INFINITE_ELEMENTS
		  // infinite elements should be skipped...
		  if (!elem->infinite())
#endif
		    { 
		      libmesh_assert (nodal_soln.size() == elem->n_nodes());
		  
		      for (unsigned int n=0; n<elem->n_nodes(); n++)
			{
			  soln[nv*(nn++) + (var + var_num)] +=
			    nodal_soln[n];
			}
		    }
		}
	    }	 
	}

      var_num += nv_sys;
    }
}



bool EquationSystems::compare (const EquationSystems& other_es, 
			       const Real threshold,
			       const bool verbose) const
{
  // safety check, whether we handle at least the same number
  // of systems
  std::vector<bool> os_result;

  if (this->n_systems() != other_es.n_systems())
    {
      if (verbose)
        {
	  std::cout << "  Fatal difference. This system handles " 
		    << this->n_systems() << " systems," << std::endl
		    << "  while the other system handles "
		    << other_es.n_systems() 
		    << " systems." << std::endl
		    << "  Aborting comparison." << std::endl;
	}
      return false;
    }
  else
    {
      // start comparing each system      
      const_system_iterator       pos = _systems.begin();
      const const_system_iterator end = _systems.end();

      for (; pos != end; ++pos)
        {
	  const std::string& sys_name = pos->first;
	  const System&  system        = *(pos->second);
      
	  // get the other system
	  const System& other_system   = other_es.get_system (sys_name);

	  os_result.push_back (system.compare (other_system, threshold, verbose));

	}

    }
  

  // sum up the results
  if (os_result.size()==0)
    return true;
  else
    {
      bool os_identical;
      unsigned int n = 0;
	  do
	    {
	      os_identical = os_result[n];
	      n++;
	    }
	  while (os_identical && n<os_result.size());
	  return os_identical;
	}
}



std::string EquationSystems::get_info () const
{
  std::ostringstream out;
  
  out << " EquationSystems\n"
      << "  n_systems()=" << this->n_systems() << '\n';

  // Print the info for the individual systems
  const_system_iterator       pos = _systems.begin();
  const const_system_iterator end = _systems.end();

  for (; pos != end; ++pos)
    out << pos->second->get_info();

  
//   // Possibly print the parameters  
//   if (!this->parameters.empty())
//     {  
//       out << "  n_parameters()=" << this->n_parameters() << '\n';
//       out << "   Parameters:\n";
      
//       for (std::map<std::string, Real>::const_iterator
// 	     param = _parameters.begin(); param != _parameters.end();
// 	   ++param)
// 	out << "    "
// 	    << "\""
// 	    << param->first
// 	    << "\""
// 	    << "="
// 	    << param->second
// 	    << '\n';
//     }
  
  return out.str();
}



void EquationSystems::print_info (std::ostream& os) const
{
  os << this->get_info()
     << std::endl;
}



std::ostream& operator << (std::ostream& os, const EquationSystems& es)
{
  es.print_info(os);
  return os;
}



unsigned int EquationSystems::n_vars () const
{
  unsigned int tot=0;
  
  const_system_iterator       pos = _systems.begin();
  const const_system_iterator end = _systems.end();

  for (; pos != end; ++pos)
    tot += pos->second->n_vars();

  return tot;
}



unsigned int EquationSystems::n_dofs () const
{
  unsigned int tot=0;

  const_system_iterator       pos = _systems.begin();
  const const_system_iterator end = _systems.end();

  for (; pos != end; ++pos)
    tot += pos->second->n_dofs();

  return tot;      
}




unsigned int EquationSystems::n_active_dofs () const
{
  unsigned int tot=0;

  const_system_iterator       pos = _systems.begin();
  const const_system_iterator end = _systems.end();

  for (; pos != end; ++pos)
    tot += pos->second->n_active_dofs();

  return tot;      
}


void EquationSystems::_add_system_to_nodes_and_elems()
{
  // All the nodes
  MeshBase::node_iterator       node_it  = _mesh.nodes_begin();
  const MeshBase::node_iterator node_end = _mesh.nodes_end();
 
  for ( ; node_it != node_end; ++node_it)
    (*node_it)->add_system();
 
  // All the elements
  MeshBase::element_iterator       elem_it  = _mesh.elements_begin();
  const MeshBase::element_iterator elem_end = _mesh.elements_end();
 
  for ( ; elem_it != elem_end; ++elem_it)
    (*elem_it)->add_system();
}