SparsityPattern::Build Class Reference

#include <dof_map.h>

List of all members.

Public Member Functions
	Build (const MeshBase &mesh_in, const DofMap &dof_map_in, const CouplingMatrix *dof_coupling_in, const bool implicit_neighbor_dofs_in, const bool need_full_sparsity_pattern_in)
	Build (Build &other, Threads::split)
void	operator() (const ConstElemRange &range)
void	join (const Build &other)
Public Attributes
SparsityPattern::Graph	sparsity_pattern
std::vector< unsigned int >	n_nz
std::vector< unsigned int >	n_oz
Private Attributes
const MeshBase &	mesh
const DofMap &	dof_map
const CouplingMatrix *	dof_coupling
const bool	implicit_neighbor_dofs
const bool	need_full_sparsity_pattern

---

Detailed Description

This helper class can be called on multiple threads to compute the sparsity pattern (or graph) of the sparse matrix resulting from the discretization. This pattern may be used directly by a particular sparse matrix format (e.g. LaspackMatrix) or indirectly (e.g. PetscMatrix). In the latter case the number of nonzeros per row of the matrix is needed for efficient preallocation. In this case it suffices to provide estimate (but bounding) values, and in this case the threaded method can take some short-cuts for efficiency.

Definition at line 100 of file dof_map.h.

---

Constructor & Destructor Documentation

SparsityPattern::Build::Build	(	const MeshBase &	mesh_in,
		const DofMap &	dof_map_in,
		const CouplingMatrix *	dof_coupling_in,
		const bool	implicit_neighbor_dofs_in,
		const bool	need_full_sparsity_pattern_in
	)			[inline]

Definition at line 115 of file dof_map.h.

                                                     :
      mesh(mesh_in),
      dof_map(dof_map_in),
      dof_coupling(dof_coupling_in),
      implicit_neighbor_dofs(implicit_neighbor_dofs_in),
      need_full_sparsity_pattern(need_full_sparsity_pattern_in)
    {}

SparsityPattern::Build::Build	(	Build &	other,
		Threads::split
	)			[inline]

Definition at line 127 of file dof_map.h.

                                       :
      mesh(other.mesh),
      dof_map(other.dof_map),
      dof_coupling(other.dof_coupling),
      implicit_neighbor_dofs(other.implicit_neighbor_dofs),
      need_full_sparsity_pattern(other.need_full_sparsity_pattern)
    {}

---

Member Function Documentation

void SparsityPattern::Build::join

(

const Build &

other

)

Definition at line 2131 of file dof_map.C.

References dof_map, DofMap::end_dof(), DofMap::first_dof(), mesh, std::min(), DofMap::n_dofs(), DofMap::n_dofs_on_processor(), n_nz, n_oz, need_full_sparsity_pattern, MeshBase::processor_id(), and sparsity_pattern.

{
  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]++;
        }
    }
}

void SparsityPattern::Build::operator()

(

const ConstElemRange &

range

)

Definition at line 1814 of file dof_map.C.

References Elem::active_family_tree_by_neighbor(), StoredRange< iterator_type, object_type >::begin(), dof_coupling, DofMap::dof_indices(), dof_map, CouplingMatrix::empty(), StoredRange< iterator_type, object_type >::end(), DofMap::end_dof(), DofMap::find_connected_dofs(), DofMap::first_dof(), implicit_neighbor_dofs, mesh, DofMap::n_dofs_on_processor(), n_nz, n_oz, Elem::n_sides(), DofMap::n_variables(), Elem::neighbor(), MeshBase::processor_id(), CouplingMatrix::size(), SparsityPattern::sort_row(), and sparsity_pattern.

{
  // 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(LIBMESH_ENABLE_AMR) || defined(LIBMESH_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 LIBMESH_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(LIBMESH_ENABLE_AMR) || defined(LIBMESH_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(LIBMESH_ENABLE_AMR) || defined(LIBMESH_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(LIBMESH_ENABLE_AMR) || defined(LIBMESH_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]++;
    }
}

---

Member Data Documentation

const CouplingMatrix* SparsityPattern::Build::dof_coupling [private]

Definition at line 105 of file dof_map.h.

Referenced by operator()().

const DofMap& SparsityPattern::Build::dof_map [private]

Definition at line 104 of file dof_map.h.

Referenced by join(), and operator()().

const bool SparsityPattern::Build::implicit_neighbor_dofs [private]

Definition at line 106 of file dof_map.h.

Referenced by operator()().

const MeshBase& SparsityPattern::Build::mesh [private]

Definition at line 103 of file dof_map.h.

Referenced by join(), and operator()().

std::vector<unsigned int> SparsityPattern::Build::n_nz

Definition at line 112 of file dof_map.h.

Referenced by DofMap::compute_sparsity(), join(), and operator()().

std::vector<unsigned int> SparsityPattern::Build::n_oz

Definition at line 113 of file dof_map.h.

Referenced by DofMap::compute_sparsity(), join(), and operator()().

const bool SparsityPattern::Build::need_full_sparsity_pattern [private]

Definition at line 107 of file dof_map.h.

Referenced by join().

SparsityPattern::Graph SparsityPattern::Build::sparsity_pattern

Definition at line 111 of file dof_map.h.

Referenced by DofMap::compute_sparsity(), join(), and operator()().

---

The documentation for this class was generated from the following files:

• dof_map.h
• dof_map.C