dof_map.c

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

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



void DofMap::dof_indices (const Elem* const elem,
			  std::vector<unsigned int>& di,
			  const unsigned int vn) const
{
  START_LOG("dof_indices()", "DofMap");
  
  libmesh_assert (elem != NULL);
  
  const unsigned int n_nodes = elem->n_nodes();
  const ElemType type        = elem->type();
  const unsigned int sys_num = this->sys_number();
  const unsigned int n_vars  = this->n_variables();
  const unsigned int dim     = elem->dim();
  
  // Clear the DOF indices vector
  di.clear();
  
#ifdef DEBUG
  // Check that sizes match in DEBUG mode
  unsigned int tot_size = 0;
#endif

  
  // Get the dof numbers
  for (unsigned int v=0; v<n_vars; v++)
    if ((v == vn) || (vn == libMesh::invalid_uint))
      { // Do this for all the variables if one was not specified
	// or just for the specified variable

        // Increase the polynomial order on p refined elements
	FEType fe_type = this->variable_type(v);
	fe_type.order = static_cast<Order>(fe_type.order +
                                           elem->p_level());

	const bool extra_hanging_dofs =
	  FEInterface::extra_hanging_dofs(fe_type);

#ifdef DEBUG
	tot_size += FEInterface::n_dofs(dim,
					fe_type,
					type);
#endif
	
	// Get the node-based DOF numbers
	for (unsigned int n=0; n<n_nodes; n++)
	  {	    
	    const Node* node      = elem->get_node(n);
	    
	    // There is a potential problem with h refinement.  Imagine a
	    // quad9 that has a linear FE on it.  Then, on the hanging side,
	    // it can falsely identify a DOF at the mid-edge node. This is why
	    // we call FEInterface instead of node->n_comp() directly.
	    const unsigned int nc = FEInterface::n_dofs_at_node (dim,
								 fe_type,
								 type,
								 n);

	    // If this is a non-vertex on a hanging node with extra
	    // degrees of freedom, we use the non-vertex dofs (which
	    // come in reverse order starting from the end, to
            // simplify p refinement)
	    if (extra_hanging_dofs && !elem->is_vertex(n))
              {
		const int dof_offset = node->n_comp(sys_num,v) - nc;

	        // We should never have fewer dofs than necessary on a
	        // node unless we're getting indices on a parent element,
	        // and we should never need the indices on such a node
	        if (dof_offset < 0)
	          {
		    libmesh_assert(!elem->active());
		    di.resize(di.size() + nc, DofObject::invalid_id);
	          }
                else
	          for (int i=node->n_comp(sys_num,v)-1; i>=dof_offset; i--)
	            {
		      libmesh_assert (node->dof_number(sys_num,v,i) !=
			      DofObject::invalid_id);
		      di.push_back(node->dof_number(sys_num,v,i));
	            }
              }
	    // If this is a vertex or an element without extra hanging
            // dofs, our dofs come in forward order coming from the
            // beginning
            else
	      for (unsigned int i=0; i<nc; i++)
	        {
		  libmesh_assert (node->dof_number(sys_num,v,i) !=
			  DofObject::invalid_id);
		  di.push_back(node->dof_number(sys_num,v,i));
	        }
	  }
	
	// If there are any element-based DOF numbers, get them
	const unsigned int nc = FEInterface::n_dofs_per_elem(dim,
							     fe_type,
							     type);
	// We should never have fewer dofs than necessary on an
	// element unless we're getting indices on a parent element,
	// and we should never need those indices
	if (nc != 0)
	  {
	    if (elem->n_systems() > sys_num &&
		nc <= elem->n_comp(sys_num,v))
	      {
		for (unsigned int i=0; i<nc; i++)
		  {
		    libmesh_assert (elem->dof_number(sys_num,v,i) !=
			    DofObject::invalid_id);
		    
		    di.push_back(elem->dof_number(sys_num,v,i));
		  }
	      }
	    else
	      {
		libmesh_assert(!elem->active() || fe_type.family == LAGRANGE);
		di.resize(di.size() + nc, DofObject::invalid_id);
	      }
	  }
      }

#ifdef DEBUG
  libmesh_assert (tot_size == di.size());
#endif
  
  STOP_LOG("dof_indices()", "DofMap");  
}
 


#ifdef ENABLE_AMR

void DofMap::old_dof_indices (const Elem* const elem,
			      std::vector<unsigned int>& di,
			      const unsigned int vn) const
{
  START_LOG("old_dof_indices()", "DofMap");
  
  libmesh_assert (elem != NULL);
  libmesh_assert (elem->old_dof_object != NULL);
	    

  const unsigned int n_nodes = elem->n_nodes();
  const ElemType type        = elem->type();
  const unsigned int sys_num = this->sys_number();
  const unsigned int n_vars  = this->n_variables();
  const unsigned int dim     = elem->dim();
  
  // Clear the DOF indices vector.
  di.clear();

#ifdef DEBUG
  // Check that sizes match
  unsigned int tot_size = 0;
#endif
  
  // Get the dof numbers
  for (unsigned int v=0; v<n_vars; v++)
    if ((v == vn) || (vn == libMesh::invalid_uint))
      { // Do this for all the variables if one was not specified
	// or just for the specified variable
	
        // Increase the polynomial order on p refined elements,
        // but make sure you get the right polynomial order for
        // the OLD degrees of freedom
        int p_adjustment = 0;
        if (elem->p_refinement_flag() == Elem::JUST_REFINED)
          {
            libmesh_assert (elem->p_level() > 0);
            p_adjustment = -1;
          }
        else if (elem->p_refinement_flag() == Elem::JUST_COARSENED)
          {
            p_adjustment = 1;
          }
	FEType fe_type = this->variable_type(v);
	fe_type.order = static_cast<Order>(fe_type.order +
                                           elem->p_level() +
                                           p_adjustment);

	const bool extra_hanging_dofs =
	  FEInterface::extra_hanging_dofs(fe_type);

#ifdef DEBUG
	tot_size += FEInterface::n_dofs(dim,
					fe_type,
					type);
#endif
	
	// Get the node-based DOF numbers
	for (unsigned int n=0; n<n_nodes; n++)
	  {
	    const Node* node      = elem->get_node(n);
	    
	    // There is a potential problem with h refinement.  Imagine a
	    // quad9 that has a linear FE on it.  Then, on the hanging side,
	    // it can falsely identify a DOF at the mid-edge node. This is why
	    // we call FEInterface instead of node->n_comp() directly.
	    const unsigned int nc = FEInterface::n_dofs_at_node (dim,
								 fe_type,
								 type,
								 n);
	    libmesh_assert (node->old_dof_object != NULL);
	    
	    // If this is a non-vertex on a hanging node with extra
	    // degrees of freedom, we use the non-vertex dofs (which
	    // come in reverse order starting from the end, to
            // simplify p refinement)
	    if (extra_hanging_dofs && !elem->is_vertex(n))
              {
	        const int dof_offset = 
                  node->old_dof_object->n_comp(sys_num,v) - nc;
	    
	        // We should never have fewer dofs than necessary on a
	        // node unless we're getting indices on a parent element
                // or a just-coarsened element
	        if (dof_offset < 0)
	          {
		    libmesh_assert(!elem->active() || elem->refinement_flag() ==
                           Elem::JUST_COARSENED);
		    di.resize(di.size() + nc, DofObject::invalid_id);
	          }
	        else
	          for (int i=node->old_dof_object->n_comp(sys_num,v)-1;
                       i>=dof_offset; i--)
	            {
		      libmesh_assert (node->old_dof_object->dof_number(sys_num,v,i) !=
			      DofObject::invalid_id);
		      di.push_back(node->old_dof_object->dof_number(sys_num,v,i));
	            }
              }
	    // If this is a vertex or an element without extra hanging
            // dofs, our dofs come in forward order coming from the
            // beginning
            else
	      for (unsigned int i=0; i<nc; i++)
	        {
		  libmesh_assert (node->old_dof_object->dof_number(sys_num,v,i) !=
			  DofObject::invalid_id);
		  di.push_back(node->old_dof_object->dof_number(sys_num,v,i));
	        }
	  }
	
	// If there are any element-based DOF numbers, get them
	const unsigned int nc = FEInterface::n_dofs_per_elem(dim,
							     fe_type,
							     type);

	// We should never have fewer dofs than necessary on an
	// element unless we're getting indices on a parent element
	// or a just-coarsened element
	if (nc != 0)
	  {
	    if (elem->old_dof_object->n_systems() > sys_num &&
		nc <= elem->old_dof_object->n_comp(sys_num,v))
	      {
		libmesh_assert (elem->old_dof_object != NULL);
		
		for (unsigned int i=0; i<nc; i++)
		  {
		    libmesh_assert (elem->old_dof_object->dof_number(sys_num,v,i) !=
			    DofObject::invalid_id);
		    
		    di.push_back(elem->old_dof_object->dof_number(sys_num,v,i));
		  }
	      }
	    else
	      {
		libmesh_assert(!elem->active() || fe_type.family == LAGRANGE ||
		       elem->refinement_flag() == Elem::JUST_COARSENED);
		di.resize(di.size() + nc, DofObject::invalid_id);
	      }
	  }
      }
  
#ifdef DEBUG
  libmesh_assert (tot_size == di.size());
#endif
  
  STOP_LOG("old_dof_indices()", "DofMap");  
}



void DofMap::augment_send_list_for_projection(const MeshBase& mesh)
{
  // Loop over the active local elements in the mesh.
  // If the element was just refined and its parent lives
  // on a different processor then we need to augment the
  // _send_list with the parent's DOF indices so that we
  // can access the parent's data for computing solution
  // projections, etc...

  // The DOF indices for the parent
  std::vector<unsigned int> di;

  // Flag telling us if we need to re-sort the send_list.
  // Note we won't need to re-sort unless a child with a
  // parent belonging to a different processor is found
  bool needs_sorting = false;
  
  
  MeshBase::const_element_iterator       elem_it  = mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator elem_end = mesh.active_local_elements_end(); 

  for ( ; elem_it != elem_end; ++elem_it)
    {
      const Elem* const elem   = *elem_it;

      // We only need to consider the children that
      // were just refined
      if (elem->refinement_flag() != Elem::JUST_REFINED) continue;
      
      const Elem* const parent = elem->parent();

      // If the parent lives on another processor
      // than the child
      if (parent != NULL)
	if (parent->processor_id() != elem->processor_id())
	  {
	    // Get the DOF indices for the parent
	    this->dof_indices (parent, di);

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

	    // We will need to re-sort the send list
	    needs_sorting = true;	    
	  }
    }

  // The send-list might need to be sorted again.
  if (needs_sorting) this->sort_send_list ();
}

#endif // ENABLE_AMR


#if defined(ENABLE_AMR) || defined(ENABLE_PERIODIC)

void DofMap::find_connected_dofs (std::vector<unsigned int>& elem_dofs) const
{

  typedef std::set<unsigned int> RCSet;

  // First insert the DOFS we already depend on into the set.  
  RCSet dof_set (elem_dofs.begin(), elem_dofs.end());

  bool done = true;
    
  // Next insert any dofs those might be constrained in terms
  // of.  Note that in this case we may not be done:  Those may
  // in turn depend on others.  So, we need to repeat this process
  // in that case until the system depends only on unconstrained
  // degrees of freedom.
  for (unsigned int i=0; i<elem_dofs.size(); i++)
    if (this->is_constrained_dof(elem_dofs[i]))
      {
	// If the DOF is constrained
	DofConstraints::const_iterator
	  pos = _dof_constraints.find(elem_dofs[i]);
	
	libmesh_assert (pos != _dof_constraints.end());
	
	const DofConstraintRow& constraint_row = pos->second;
	
// adaptive p refinement currently gives us lots of empty constraint
// rows - we should optimize those DoFs away in the future.  [RHS]
//	libmesh_assert (!constraint_row.empty());
	
	DofConstraintRow::const_iterator it     = constraint_row.begin();
	DofConstraintRow::const_iterator it_end = constraint_row.end();
	
	
	// Add the DOFs this dof is constrained in terms of.
	// note that these dofs might also be constrained, so
	// we will need to call this function recursively.
	for ( ; it != it_end; ++it)
	  if (!dof_set.count (it->first))
	    {
	      dof_set.insert (it->first);
	      done = false;
	    }
	}
  
  
  // If not done then we need to do more work
  // (obviously :-) )!
  if (!done)
    {
      // Fill the vector with the contents of the set
      elem_dofs.clear();