dof_map.c

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

          unsigned int n_comp =
            obj->n_comp(this->sys_number(), v);
          unsigned int first_dof = n_comp ?
            obj->dof_number(this->sys_number(), v, 0) : 0;
          libmesh_assert(first_dof != DofObject::invalid_id);
        }
    }
#endif // DEBUG

  if (!build_send_list)
    {
      _send_list.clear();
      _send_list.reserve(send_list_size);
    }
}



void DofMap::distribute_local_dofs_var_major(unsigned int &next_free_dof,
                                             MeshBase& mesh,
                                             const bool build_send_list)
{
  const unsigned int sys_num = this->sys_number();
  const unsigned int n_vars  = this->n_variables();

  // Number of elements to add to send_list
  unsigned int send_list_size = 0;
  
  //-------------------------------------------------------------------------
  // First count and assign temporary numbers to local dofs
  for (unsigned var=0; var<n_vars; var++)
    {
      MeshBase::element_iterator       elem_it  = mesh.active_local_elements_begin();
      const MeshBase::element_iterator elem_end = mesh.active_local_elements_end();

      for ( ; elem_it != elem_end; ++elem_it)
        {
          // Only number dofs connected to active
          // elements on this processor.
          Elem* elem                 = *elem_it;
          const unsigned int n_nodes = elem->n_nodes();
          
          // First number the nodal DOFS
          for (unsigned int n=0; n<n_nodes; n++)
            {
              Node* node = elem->get_node(n);
              
              // assign dof numbers (all at once) if this is
              // our node and if they aren't already there
              if ((node->n_comp(sys_num,var) > 0) &&
                  (node->processor_id() == libMesh::processor_id()) &&
                  (node->dof_number(sys_num,var,0) ==
                   DofObject::invalid_id))
                {
                  node->set_dof_number(sys_num,
                                       var,
                                       0,
                                       next_free_dof);
                  next_free_dof += node->n_comp(sys_num,var);

		  // If we've been requested to, add local dofs to the
		  // _send_list
                  if (build_send_list)
                    for (unsigned int index=0; index<node->n_comp(sys_num,var);
                         index++)
                      _send_list.push_back(node->dof_number(sys_num,
                                                            var,
                                                            index));
                  else
                    send_list_size += node->n_comp(sys_num,var);
                }
            }
                  
          // Now number the element DOFS
          if (elem->n_comp(sys_num,var) > 0)
            {
              libmesh_assert (elem->dof_number(sys_num,var,0) ==
                      DofObject::invalid_id);

              elem->set_dof_number(sys_num,
                                   var,
                                   0,
                                   next_free_dof);

              next_free_dof += elem->n_comp(sys_num,var);
              
	      // If we've been requested to, add local dofs to the
	      // _send_list
              if (build_send_list)
                for (unsigned int index=0; index<elem->n_comp(sys_num,var);
                     index++)
                  _send_list.push_back(elem->dof_number(sys_num,
                                                      var,
                                                      index));
              else
                send_list_size += elem->n_comp(sys_num,var);
            }
        }
    }

#ifdef DEBUG
// Make sure we didn't miss any nodes
  MeshTools::libmesh_assert_valid_node_procids(mesh);

  MeshBase::node_iterator       node_it  = mesh.local_nodes_begin();
  const MeshBase::node_iterator node_end = mesh.local_nodes_end();
  for (; node_it != node_end; ++node_it)
    {
      Node *obj = *node_it;
      libmesh_assert(obj);
      unsigned int n_variables = obj->n_vars(this->sys_number());
      for (unsigned int v=0; v != n_variables; ++v)
        {
          unsigned int n_comp =
            obj->n_comp(this->sys_number(), v);
          unsigned int first_dof = n_comp ?
            obj->dof_number(this->sys_number(), v, 0) : 0;
          libmesh_assert(first_dof != DofObject::invalid_id);
        }
    }
#endif // DEBUG

  if (!build_send_list)
    {
      _send_list.clear();
      _send_list.reserve(send_list_size);
    }
}



void DofMap::add_neighbors_to_send_list(MeshBase& mesh)
{
  START_LOG("add_neighbors_to_send_list()", "DofMap");
  
  //-------------------------------------------------------------------------
  // The _send_list now only contains entries from the local processor.
  // We need to add the DOFs from elements that live on neighboring processors
  // that are neighbors of the elements on the local processor
  //-------------------------------------------------------------------------

  MeshBase::const_element_iterator       local_elem_it
    = mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator local_elem_end
    = mesh.active_local_elements_end(); 

  std::vector<bool> node_on_processor(mesh.max_node_id(), false);
  std::vector<unsigned int> di;
  std::vector<const Elem *> family;

  // Loop over the active local elements, adding all active elements
  // that neighbor an active local element to the send list.
  for ( ; local_elem_it != local_elem_end; ++local_elem_it)
    {
      const Elem* elem = *local_elem_it;

      // Flag all the nodes of active local elements as seen
      for (unsigned int n=0; n!=elem->n_nodes(); n++)
        node_on_processor[elem->node(n)] = true;

      // Loop over the neighbors of those elements
      for (unsigned int s=0; s<elem->n_neighbors(); s++)
	if (elem->neighbor(s) != NULL)
	  {
	    family.clear();
	    
            // Find all the active elements that neighbor elem
#ifdef ENABLE_AMR
            if (!elem->neighbor(s)->active())
              elem->neighbor(s)->active_family_tree_by_neighbor(family, elem);
            else
#endif
              family.push_back(elem->neighbor(s));

            for (unsigned int i=0; i!=family.size(); ++i)
	    // If the neighbor lives on a different processor
	    if (family[i]->processor_id() != libMesh::processor_id())
	      {
	        // Get the DOF indices for this neighboring element
	        this->dof_indices (family[i], di);

	        // Insert the DOF indices into the send list
	        _send_list.insert (_send_list.end(),
				   di.begin(), di.end());
	      }
	  }
    }

  // Now loop over all non_local active elements and add any missing
  // nodal-only neighbors
  MeshBase::const_element_iterator       elem_it
    = mesh.active_elements_begin();
  const MeshBase::const_element_iterator elem_end
    = mesh.active_elements_end();
  
  for ( ; elem_it != elem_end; ++elem_it)
    {
      const Elem* elem = *elem_it;

      // If this is one of our elements, we've already added it
      if (elem->processor_id() == libMesh::processor_id())
        continue;

      // Do we need to add the element DOFs?
      bool add_elem_dofs = false;
      
      // Check all the nodes of the element to see if it
      // shares a node with us
      for (unsigned int n=0; n!=elem->n_nodes(); n++)
        if (node_on_processor[elem->node(n)])
	  add_elem_dofs = true;

      // Add the element degrees of freedom if it shares at
      // least one node.
      if (add_elem_dofs)
	{
	  // Get the DOF indices for this neighboring element
	  this->dof_indices (elem, di);

	  // Insert the DOF indices into the send list
	  _send_list.insert (_send_list.end(),
			       di.begin(), di.end());
	}
    }
  
  STOP_LOG("add_neighbors_to_send_list()", "DofMap");
}



void DofMap::sort_send_list ()
{
  START_LOG("sort_send_list()", "DofMap");
  
  // First sort the send list.  After this
  // duplicated elements will be adjacent in the
  // vector
  std::sort(_send_list.begin(), _send_list.end());

  // Now use std::unique to remove duplicate entries  
  std::vector<unsigned int>::iterator new_end =
    std::unique (_send_list.begin(), _send_list.end());

  // Remove the end of the send_list.  Use the "swap trick"
  // from Effective STL
  std::vector<unsigned int> (_send_list.begin(), new_end).swap (_send_list);
  
  STOP_LOG("sort_send_list()", "DofMap");
}



void DofMap::compute_sparsity(const MeshBase& mesh)
{
  libmesh_assert (mesh.is_prepared());
  libmesh_assert (this->n_variables());

  START_LOG("compute_sparsity()", "DofMap");

  // Compute the sparsity structure of the global matrix.  This can be
  // fed into a PetscMatrix to allocate exacly the number of nonzeros
  // necessary to store the matrix.  This algorithm should be linear
  // in the (# of elements)*(# nodes per element)

  bool implicit_neighbor_dofs = 
    libMesh::on_command_line ("--implicit_neighbor_dofs");

  // look at all the variables in this system.  If any are discontinuous then
  // force implicit_neighbor_dofs=true
  {
    bool all_discontinuous_dofs = true;
    
    for (unsigned int var=0; var<this->n_variables(); var++)
      if (FEBase::build (mesh.mesh_dimension(),
			 this->variable_type(var))->get_continuity() !=  DISCONTINUOUS)
	all_discontinuous_dofs = false;

    if (all_discontinuous_dofs)
      implicit_neighbor_dofs = true;
  }
  
  // We can be more efficient in the threaded sparsity pattern assembly
  // if we don't need the exact pattern.  For some sparse matrix formats
  // a good upper bound will suffice.
  bool need_full_sparsity_pattern=false;
  std::vector<SparseMatrix<Number>* >::iterator
    pos = _matrices.begin(),
    end = _matrices.end();
      
  for (; pos != end; ++pos)
    if ((*pos)->need_full_sparsity_pattern())
      need_full_sparsity_pattern = true;
  
  // We can compute the sparsity pattern in parallel on multiple
  // threads.  The goal is for each thread to compute the full sparsity
  // pattern for a subset of elements.  These sparsity patterns can
  // be efficiently merged in the SparsityPattern::Build::join()
  // method, especially if there is not too much overlap between them.
  // Even better, if the full sparsity pattern is not needed then
  // the number of nonzeros per row can be estimated from the
  // sparsity patterns created on each thread.
  SparsityPattern::Build sp (mesh,
			     *this,
			     _dof_coupling,
			     implicit_neighbor_dofs,
			     need_full_sparsity_pattern);
  
  Threads::parallel_reduce (ConstElemRange (mesh.active_elements_begin(),
					    mesh.active_elements_end()), sp);

#ifndef NDEBUG
  // Avoid declaring these variables unless asserts are enabled.
  const unsigned int proc_id        = mesh.processor_id();
  const unsigned int n_dofs_on_proc = this->n_dofs_on_processor(proc_id);
#endif
  libmesh_assert (sp.sparsity_pattern.size() == n_dofs_on_proc);

  // steal the n_nz and n_oz arrays from sp -- it won't need them any more,
  // and this is more efficient than copying them.
  _n_nz.swap(sp.n_nz);
  _n_oz.swap(sp.n_oz);
  
  STOP_LOG("compute_sparsity()", "DofMap");
  
  // We are done with the sparsity_pattern.  However, quite a
  // lot has gone into computing it.  It is possible that some
  // \p SparseMatrix implementations want to see it.  Let them
  // see it before we throw it away.
  pos = _matrices.begin();
  end = _matrices.end();
      
  for (; pos != end; ++pos)
    (*pos)->update_sparsity_pattern (sp.sparsity_pattern);     
}



void DofMap::extract_local_vector (const NumericVector<Number>& Ug,
				   const std::vector<unsigned int>& dof_indices,
				   DenseVectorBase<Number>& Ue) const
{
#ifdef ENABLE_AMR

  // Trivial mapping
  libmesh_assert (dof_indices.size() == Ue.size());
  bool has_constrained_dofs = false;

  for (unsigned int il=0; il<dof_indices.size(); il++)
    {
      const unsigned int ig = dof_indices[il];

      if (this->is_constrained_dof (ig)) has_constrained_dofs = true;
      
      libmesh_assert ((il >= Ug.first_local_index()) &&
	      (il <  Ug.last_local_index()));

      Ue.el(il) = Ug(ig);
    }

  // If the element has any constrained DOFs then we need
  // to account for them in the mapping.  This will handle
  // the case that the input vector is not constrained.
  if (has_constrained_dofs)
    {
      // Copy the input DOF indices.
      std::vector<unsigned int> constrained_dof_indices(dof_indices);

      DenseMatrix<Number> C;

      this->build_constraint_matrix (C, constrained_dof_indices);

      libmesh_assert (dof_indices.size()             == C.m());
      libmesh_assert (constrained_dof_indices.size() == C.n());

      // zero-out Ue
      Ue.zero();

      // compute Ue = C Ug, with proper mapping.
      for (unsigned int i=0; i<dof_indices.size(); i++)
	for (unsigned int j=0; j<constrained_dof_indices.size(); j++)
	  {
	    const unsigned int jg = constrained_dof_indices[j];

	    libmesh_assert ((jg >= Ug.first_local_index()) &&
		    (jg <  Ug.last_local_index()));
	    
	    Ue.el(i) += C(i,j)*Ug(jg); 	    
	  }
    }  
   
#else
  
  // Trivial mapping
  libmesh_assert (dof_indices.size() == Ue.size());
  
  for (unsigned int il=0; il<dof_indices.size(); il++)
    {
      const unsigned int ig = dof_indices[il];
      
      libmesh_assert ((ig >= Ug.first_local_index()) && (ig <  Ug.last_local_index()));