dof_map.c

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

      elem_dofs.insert (elem_dofs.end(),
			dof_set.begin(), dof_set.end());
      

      // May need to do this recursively.  It is possible
      // that we just replaced a constrained DOF with another
      // constrained DOF.
      this->find_connected_dofs (elem_dofs);
      
    } // end if (!done)
}

#endif // ENABLE_AMR || ENABLE_PERIODIC



#if defined(__GNUC__) && (__GNUC__ < 4) && !defined(__INTEL_COMPILER)

void SparsityPattern::_dummy_function(void)
{
}

#endif



void SparsityPattern::Build::operator()(const ConstElemRange &range)
{
  // 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)
  const unsigned int proc_id           = mesh.processor_id();
  const unsigned int n_dofs_on_proc    = dof_map.n_dofs_on_processor(proc_id);
  const unsigned int first_dof_on_proc = dof_map.first_dof(proc_id);
  const unsigned int end_dof_on_proc   = dof_map.end_dof(proc_id);

  sparsity_pattern.resize(n_dofs_on_proc);
  
  // If the user did not explicitly specify the DOF coupling
  // then all the DOFS are coupled to each other.  Furthermore,
  // we can take a shortcut and do this more quickly here.  So
  // we use an if-test.
  if ((dof_coupling == NULL) || (dof_coupling->empty()))
    {
      std::vector<unsigned int>
	element_dofs,
	neighbor_dofs,
	dofs_to_add;

      std::vector<const Elem*> active_neighbors;

      for (ConstElemRange::const_iterator elem_it = range.begin() ; elem_it != range.end(); ++elem_it)
	{
	  const Elem* const elem = *elem_it;

	  // Get the global indices of the DOFs with support on this element
	  dof_map.dof_indices (elem, element_dofs);
#if defined(ENABLE_AMR) || defined(ENABLE_PERIODIC)
	  dof_map.find_connected_dofs (element_dofs);
#endif

	  // We can be more efficient if we sort the element DOFs
	  // into increasing order
	  std::sort(element_dofs.begin(), element_dofs.end());
	  
	  const unsigned int n_dofs_on_element = element_dofs.size();
	    
	  for (unsigned int i=0; i<n_dofs_on_element; i++)
	    {
	      const unsigned int ig = element_dofs[i];
	      
	      // Only bother if this matrix row will be stored
	      // on this processor.
	      if ((ig >= first_dof_on_proc) &&
		  (ig <  end_dof_on_proc))
		{
		  // This is what I mean
		  // libmesh_assert ((ig - first_dof_on_proc) >= 0);
		  // but do the test like this because ig and
		  // first_dof_on_proc are unsigned ints
		  libmesh_assert (ig >= first_dof_on_proc);
		  libmesh_assert ((ig - first_dof_on_proc) < sparsity_pattern.size());
		  
		  SparsityPattern::Row &row = sparsity_pattern[ig - first_dof_on_proc];

		  // If the row is empty we will add *all* the element DOFs,
		  // so just do that.
		  if (row.empty())
		    {
 		      row.insert(row.end(),
				 element_dofs.begin(),
				 element_dofs.end());
		    }
		  else
		    {		  
		      // Build a list of the DOF indices not found in the
		      // sparsity pattern
		      dofs_to_add.clear();

		      // Cache iterators.  Low will move forward, subsequent
		      // searches will be on smaller ranges
		      SparsityPattern::Row::iterator
			low  = std::lower_bound (row.begin(), row.end(), element_dofs.front()),
			high = std::upper_bound (low,         row.end(), element_dofs.back());
		  
		      for (unsigned int j=0; j<n_dofs_on_element; j++)
			{
			  const unsigned int jg = element_dofs[j];
			  
			  // See if jg is in the sorted range
			  std::pair<SparsityPattern::Row::iterator,
			            SparsityPattern::Row::iterator>
			    pos = std::equal_range (low, high, jg);
			
			  // Must add jg if it wasn't found
			  if (pos.first == pos.second)
			    dofs_to_add.push_back(jg);
			
			  // pos.first is now a valid lower bound for any
			  // remaining element DOFs. (That's why we sorted them.)
			  // Use it for the next search
			  low = pos.first;
			}

		      // Add to the sparsity pattern
		      if (!dofs_to_add.empty())
			{
			  const unsigned int old_size = row.size();
			  
			  row.insert (row.end(),
				      dofs_to_add.begin(),
				      dofs_to_add.end());
		      
			  SparsityPattern::sort_row (row.begin(), row.begin()+old_size, row.end());
			}
		    }
		  
		  // Now (possibly) add dofs from neighboring elements
		  // TODO:[BSK] optimize this like above!
		  if (implicit_neighbor_dofs)
		    for (unsigned int s=0; s<elem->n_sides(); s++)
		      if (elem->neighbor(s) != NULL)
			{
			  const Elem* const neighbor_0 = elem->neighbor(s);
#ifdef ENABLE_AMR
			  neighbor_0->active_family_tree_by_neighbor(active_neighbors,elem);
#else
			  active_neighbors.clear();
			  active_neighbors.push_back(neighbor_0);
#endif
			  
			  for (unsigned int a=0; a != active_neighbors.size(); ++a)
			    {
			      const Elem *neighbor = active_neighbors[a];
			      
			      dof_map.dof_indices (neighbor, neighbor_dofs);
#if defined(ENABLE_AMR) || defined(ENABLE_PERIODIC)
			      dof_map.find_connected_dofs (neighbor_dofs);
#endif			      
			      const unsigned int n_dofs_on_neighbor = neighbor_dofs.size();
			      
			      for (unsigned int j=0; j<n_dofs_on_neighbor; j++)
				{
				  const unsigned int jg = neighbor_dofs[j];
				  
				  // See if jg is in the sorted range
				  std::pair<SparsityPattern::Row::iterator,
				            SparsityPattern::Row::iterator>
				    pos = std::equal_range (row.begin(), row.end(), jg);
			        
				  // Insert jg if it wasn't found
				  if (pos.first == pos.second)
				    row.insert (pos.first, jg);
				}    
			    }
			}
		}
	    }
	}	      
    } 


  
  // This is what we do in the case that the user has specified
  // explicit DOF coupling.
  else
    {
      libmesh_assert (dof_coupling != NULL);
      libmesh_assert (dof_coupling->size() ==
	      dof_map.n_variables());
      
      const unsigned int n_var = dof_map.n_variables();
      
      std::vector<unsigned int>
	element_dofs_i,
	element_dofs_j,
	dofs_to_add;


      for (ConstElemRange::const_iterator elem_it = range.begin() ; elem_it != range.end(); ++elem_it)
	for (unsigned int vi=0; vi<n_var; vi++)
	  {
	    const Elem* const elem = *elem_it;
	    
	    // Find element dofs for variable vi
	    dof_map.dof_indices (elem, element_dofs_i, vi);
#if defined(ENABLE_AMR) || defined(ENABLE_PERIODIC)
	    dof_map.find_connected_dofs (element_dofs_i);
#endif

	    // We can be more efficient if we sort the element DOFs
	    // into increasing order
	    std::sort(element_dofs_i.begin(), element_dofs_i.end());
	    const unsigned int n_dofs_on_element_i = element_dofs_i.size();
	    
	    for (unsigned int vj=0; vj<n_var; vj++)
	      if ((*dof_coupling)(vi,vj)) // If vi couples to vj
		{
		  // Find element dofs for variable vj, note that
		  // if vi==vj we already have the dofs.
		  if (vi != vj)
		    {
		      dof_map.dof_indices (elem, element_dofs_j, vj);
#if defined(ENABLE_AMR) || defined(ENABLE_PERIODIC)
		      dof_map.find_connected_dofs (element_dofs_j);
#endif

		      // We can be more efficient if we sort the element DOFs
		      // into increasing order
		      std::sort (element_dofs_j.begin(), element_dofs_j.end());
		    }
		  else
		    element_dofs_j = element_dofs_i;
		  
		  const unsigned int n_dofs_on_element_j =
		    element_dofs_j.size();
		  
		  for (unsigned int i=0; i<n_dofs_on_element_i; i++)
		    {
		      const unsigned int ig = element_dofs_i[i];
		      
		      // We are only interested in computing the sparsity pattern
		      // for the rows in [range.begin(),range.end()).  So, if this
		      // element's DOFs are not in that range just go ahead to the
		      // next one.
		      //	      
		      // Also, only bother if this matrix row will be stored
		      // on this processor.
		      if ((ig >= first_dof_on_proc) &&
			  (ig <  end_dof_on_proc))			
			{
			  // This is what I mean
			  // libmesh_assert ((ig - first_dof_on_proc) >= 0);
			  // but do the test like this because ig and
			  // first_dof_on_proc are unsigned ints
			  libmesh_assert (ig >= first_dof_on_proc);
			  libmesh_assert (ig < (sparsity_pattern.size() +
					first_dof_on_proc));
			  
			  SparsityPattern::Row &row = sparsity_pattern[ig - first_dof_on_proc];

			  // If the row is empty we will add *all* the element j DOFs,
			  // so just do that.
			  if (row.empty())
			    {
			      row.insert(row.end(),
					 element_dofs_j.begin(),
					 element_dofs_j.end());
			    }
			  else
			    {		  
			      // Build a list of the DOF indices not found in the
			      // sparsity pattern
			      dofs_to_add.clear();

			      // Cache iterators.  Low will move forward, subsequent
			      // searches will be on smaller ranges
			      SparsityPattern::Row::iterator
				low  = std::lower_bound (row.begin(), row.end(), element_dofs_j.front()),
				high = std::upper_bound (low,         row.end(), element_dofs_j.back());

			      for (unsigned int j=0; j<n_dofs_on_element_j; j++)
				{
				  const unsigned int jg = element_dofs_j[j];
				  
				  // See if jg is in the sorted range
				  std::pair<SparsityPattern::Row::iterator,
				            SparsityPattern::Row::iterator>
				    pos = std::equal_range (low, high, jg);
			      
				  // Must add jg if it wasn't found
				  if (pos.first == pos.second)
				    dofs_to_add.push_back(jg);
				
				  // pos.first is now a valid lower bound for any
				  // remaining element j DOFs. (That's why we sorted them.)
				  // Use it for the next search
				  low = pos.first;
				}
			      
			      // Add to the sparsity pattern
			      if (!dofs_to_add.empty())
				{
				  const unsigned int old_size = row.size();
				  
				  row.insert (row.end(),
					      dofs_to_add.begin(),
					      dofs_to_add.end());
		      
				  SparsityPattern::sort_row (row.begin(), row.begin()+old_size, row.end());
				}
			    }
			}
		    }
		}
	  }
    }
  
  // Now the full sparsity structure is built for all of the
  // DOFs connected to our rows of the matrix.
  n_nz.resize (n_dofs_on_proc);
  n_oz.resize (n_dofs_on_proc);

  // First zero the counters.
  std::fill(n_nz.begin(), n_nz.end(), 0);
  std::fill(n_oz.begin(), n_oz.end(), 0);
  
  for (unsigned int i=0; i<n_dofs_on_proc; i++)
    {
      // Get the row of the sparsity pattern
      const SparsityPattern::Row &row = sparsity_pattern[i];

      for (unsigned int j=0; j<row.size(); j++)
	if ((row[j] < first_dof_on_proc) || (row[j] >= end_dof_on_proc))
	  n_oz[i]++;
	else
	  n_nz[i]++;
    }
}



void SparsityPattern::Build::join (const SparsityPattern::Build &other)
{
  const unsigned int proc_id           = mesh.processor_id();
  const unsigned int n_global_dofs     = dof_map.n_dofs();
  const unsigned int n_dofs_on_proc    = dof_map.n_dofs_on_processor(proc_id);
  const unsigned int first_dof_on_proc = dof_map.first_dof(proc_id);
  const unsigned int end_dof_on_proc   = dof_map.end_dof(proc_id);

  libmesh_assert (sparsity_pattern.size() ==  other.sparsity_pattern.size());
  libmesh_assert (n_nz.size() == sparsity_pattern.size());
  libmesh_assert (n_oz.size() == sparsity_pattern.size());
  
  for (unsigned int r=0; r<n_dofs_on_proc; r++)
    {
      // incriment the number of on and off-processor nonzeros in this row
      // (note this will be an upper bound unless we need the full sparsity pattern)
      n_nz[r] += other.n_nz[r];  n_nz[r] = std::min(n_nz[r], n_dofs_on_proc);
      n_oz[r] += other.n_oz[r];  n_oz[r] = std::min(n_oz[r], n_global_dofs-n_nz[r]);

      if (need_full_sparsity_pattern)
	{
	  SparsityPattern::Row       &my_row    = sparsity_pattern[r];
	  const SparsityPattern::Row &their_row = other.sparsity_pattern[r];
      
	  // simple copy if I have no dofs
	  if (my_row.empty())
	    my_row = their_row;
	  
	  // otherwise add their DOFs to mine, resort, and re-unique the row
	  else if (!their_row.empty()) // do nothing for the trivial case where 
	    {                          // their row is empty
	      my_row.insert (my_row.end(),
			     their_row.begin(),
			     their_row.end());

	      // We cannot use SparsityPattern::sort_row() here because it expects
	      // the [begin,middle) [middle,end) to be non-overlapping.  This is not
	      // necessarily the case here, so use std::sort()
	      std::sort (my_row.begin(), my_row.end());
	 
	      my_row.erase(std::unique (my_row.begin(), my_row.end()), my_row.end());
	    }

	  // fix the number of on and off-processor nonzeros in this row
	  n_nz[r] = n_oz[r] = 0;
	  
	  for (unsigned int j=0; j<my_row.size(); j++)
	    if ((my_row[j] < first_dof_on_proc) || (my_row[j] >= end_dof_on_proc))
	      n_oz[r]++;
	    else
	      n_nz[r]++;
	}
    }
}