DofMap Class Reference

#include <dof_map.h>

Inheritance diagram for DofMap:

[legend]

List of all members.

Public Member Functions
	DofMap (const unsigned int sys_number)
	~DofMap ()
void	attach_matrix (SparseMatrix< Number > &matrix)
void	distribute_dofs (MeshBase &)
void	compute_sparsity (const MeshBase &)
void	prepare_send_list ()
const std::vector< unsigned int > &	get_send_list () const
const std::vector< unsigned int > &	get_n_nz () const
const std::vector< unsigned int > &	get_n_oz () const
void	add_variable (const System::Variable &var)
const System::Variable &	variable (const unsigned int c) const
Order	variable_order (const unsigned int c) const
const FEType &	variable_type (const unsigned int c) const
unsigned int	n_variables () const
unsigned int	n_dofs () const
unsigned int	n_SCALAR_dofs () const
unsigned int	n_local_dofs () const
unsigned int	n_dofs_on_processor (const unsigned int proc) const
unsigned int	first_dof (const unsigned int proc=libMesh::processor_id()) const
unsigned int	last_dof (const unsigned int proc=libMesh::processor_id()) const
unsigned int	end_dof (const unsigned int proc=libMesh::processor_id()) const
unsigned int	variable_first_local_dof (const unsigned int var) const
unsigned int	variable_last_local_dof (const unsigned int var) const
void	dof_indices (const Elem *const elem, std::vector< unsigned int > &di, const unsigned int vn=libMesh::invalid_uint) const
void	SCALAR_dof_indices (std::vector< unsigned int > &di, const unsigned int vn, const bool old_dofs=false) const
bool	use_coupled_neighbor_dofs (const MeshBase &mesh) const
void	extract_local_vector (const NumericVector< Number > &Ug, const std::vector< unsigned int > &dof_indices, DenseVectorBase< Number > &Ue) const
unsigned int	n_constrained_dofs () const
void	create_dof_constraints (const MeshBase &mesh)
void	allgather_recursive_constraints ()
void	process_constraints ()
void	add_constraint_row (const unsigned int dof_number, const DofConstraintRow &constraint_row, const bool forbid_constraint_overwrite=true)
bool	is_constrained_dof (const unsigned int dof) const
void	print_dof_constraints (std::ostream &os=std::cout) const
std::pair< Real, Real >	max_constraint_error (const System &system, NumericVector< Number > *v=NULL) const
void	constrain_element_matrix (DenseMatrix< Number > &matrix, std::vector< unsigned int > &elem_dofs, bool asymmetric_constraint_rows=true) const
void	constrain_element_matrix (DenseMatrix< Number > &matrix, std::vector< unsigned int > &row_dofs, std::vector< unsigned int > &col_dofs, bool asymmetric_constraint_rows=true) const
void	constrain_element_vector (DenseVector< Number > &rhs, std::vector< unsigned int > &dofs, bool asymmetric_constraint_rows=true) const
void	constrain_element_matrix_and_vector (DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< unsigned int > &elem_dofs, bool asymmetric_constraint_rows=true) const
void	constrain_element_dyad_matrix (DenseVector< Number > &v, DenseVector< Number > &w, std::vector< unsigned int > &row_dofs, bool asymmetric_constraint_rows=true) const
void	constrain_nothing (std::vector< unsigned int > &dofs) const
void	enforce_constraints_exactly (const System &system, NumericVector< Number > *v=NULL) const
void	add_periodic_boundary (const PeriodicBoundary &periodic_boundary)
bool	is_periodic_boundary (const unsigned int boundaryid) const
void	old_dof_indices (const Elem *const elem, std::vector< unsigned int > &di, const unsigned int vn=libMesh::invalid_uint) const
unsigned int	n_old_dofs () const
void	constrain_p_dofs (unsigned int var, const Elem *elem, unsigned int s, unsigned int p)
void	reinit (MeshBase &mesh)
void	clear ()
Static Public Member Functions
static std::string	get_info ()
static void	print_info ()
static unsigned int	n_objects ()
Public Attributes
CouplingMatrix *	_dof_coupling
Protected Types
typedef std::map< std::string, std::pair< unsigned int, unsigned int > >	Counts
Protected Member Functions
void	increment_constructor_count (const std::string &name)
void	increment_destructor_count (const std::string &name)
Static Protected Attributes
static Counts	_counts
static Threads::atomic < unsigned int >	_n_objects
static Threads::spin_mutex	_mutex
Private Types
typedef DofObject *(DofMap::*	dofobject_accessor )(MeshBase &mesh, unsigned int i) const
Private Member Functions
unsigned int	sys_number () const
void	invalidate_dofs (MeshBase &mesh) const
DofObject *	node_ptr (MeshBase &mesh, unsigned int i) const
DofObject *	elem_ptr (MeshBase &mesh, unsigned int i) const
template<typename iterator_type >
void	set_nonlocal_dof_objects (iterator_type objects_begin, iterator_type objects_end, MeshBase &mesh, dofobject_accessor objects)
void	distribute_local_dofs_var_major (unsigned int &next_free_dof, MeshBase &mesh)
void	distribute_local_dofs_node_major (unsigned int &next_free_dof, MeshBase &mesh)
void	add_neighbors_to_send_list (MeshBase &mesh)
void	build_constraint_matrix (DenseMatrix< Number > &C, std::vector< unsigned int > &elem_dofs, const bool called_recursively=false) const
void	find_connected_dofs (std::vector< unsigned int > &elem_dofs) const
void	add_constraints_to_send_list ()
Private Attributes
std::vector< System::Variable >	_variables
const unsigned int	_sys_number
std::vector< SparseMatrix < Number > * >	_matrices
std::vector< unsigned int >	_first_df
std::vector< unsigned int >	_end_df
std::vector< unsigned int >	_var_first_local_df
std::vector< unsigned int >	_send_list
std::vector< unsigned int >	_n_nz
std::vector< unsigned int >	_n_oz
unsigned int	_n_dfs
unsigned int	_n_SCALAR_dofs
unsigned int	_n_old_dfs
DofConstraints	_dof_constraints
PeriodicBoundaries	_periodic_boundaries
Friends
class	SparsityPattern::Build

---

Detailed Description

This class handles the numbering of degrees of freedom on a mesh. For systems of equations the class supports a fixed number of variables. The degrees of freedom are numbered such that sequential, contiguous blocks correspond to distinct subdomains. This is so that the resulting data structures will work well with parallel linear algebra packages.

Author:

Benjamin S. Kirk, 2002-2007

Definition at line 241 of file dof_map.h.

---

Member Typedef Documentation

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > ReferenceCounter::Counts [protected, inherited]

Data structure to log the information. The log is identified by the class name.

Definition at line 105 of file reference_counter.h.

typedef DofObject*(DofMap::* DofMap::dofobject_accessor)(MeshBase &mesh, unsigned int i) const [private]

A member function type like node_ptr or elem_ptr

---

Constructor & Destructor Documentation

DofMap::DofMap

(

const unsigned int

sys_number

)

[inline]

Constructor. Requires the number of the system for which we will be numbering degrees of freedom.

Definition at line 878 of file dof_map.h.

References _matrices.

                                        :
  _dof_coupling(NULL),
  _sys_number(number),
//  _matrix(NULL),
  _n_dfs(0),
  _n_SCALAR_dofs(0)
#ifdef LIBMESH_ENABLE_AMR
  , _n_old_dfs(0)
#endif
{
  _matrices.clear();
}

DofMap::~DofMap

(

)

Destructor.

Definition at line 49 of file dof_map.C.

References clear().

{
  this->clear();
}

---

Member Function Documentation

void DofMap::add_constraint_row	(	const unsigned int	dof_number,
		const DofConstraintRow &	constraint_row,
		const bool	forbid_constraint_overwrite = true
	)

Adds a copy of the user-defined row to the constraint matrix. By default, produces an error if the DOF was already constrained.

Definition at line 181 of file dof_map_constraints.C.

References _dof_constraints, and is_constrained_dof().

{
  // Optionally allow the user to overwrite constraints.  Defaults to false.
  if (forbid_constraint_overwrite)
    if (this->is_constrained_dof(dof_number))
      {
        std::cerr << "ERROR: DOF " << dof_number << " was already constrained!"
                  << std::endl;
        libmesh_error();
      }
  
  std::pair<unsigned int, DofConstraintRow> kv(dof_number, constraint_row);

  _dof_constraints.insert(kv);
}

void DofMap::add_constraints_to_send_list

(

)

[private]

Adds entries to the _send_list vector corresponding to DoFs which are dependencies for constraint equations on the current processor.

Definition at line 1192 of file dof_map_constraints.C.

References _dof_constraints, _send_list, end_dof(), first_dof(), std::max(), and libMesh::n_processors().

Referenced by process_constraints().

{
  // This function must be run on all processors at once
  parallel_only();

  // Return immediately if there's nothing to gather
  if (libMesh::n_processors() == 1)
    return;

  // We might get to return immediately if none of the processors
  // found any constraints
  unsigned int has_constraints = !_dof_constraints.empty();
  Parallel::max(has_constraints);
  if (!has_constraints)
    return;

  for (DofConstraints::iterator i = _dof_constraints.begin();
         i != _dof_constraints.end(); ++i)
    {
      unsigned int constrained_dof = i->first;

      // We only need the dependencies of our own constrained dofs
      if (constrained_dof < this->first_dof() || 
          constrained_dof >= this->end_dof())
        continue;

      DofConstraintRow& constraint_row = i->second;
      for (DofConstraintRow::const_iterator
           j=constraint_row.begin(); j != constraint_row.end();
           ++j)
        {
          unsigned int constraint_dependency = j->first;

          // No point in adding one of our own dofs to the send_list
          if (constraint_dependency >= this->first_dof() && 
              constraint_dependency < this->end_dof())
            continue;

          _send_list.push_back(constraint_dependency);
        }
    }
}

void DofMap::add_neighbors_to_send_list

(

MeshBase &

mesh

)

[private]

Adds entries to the _send_list vector corresponding to DoFs on elements neighboring the current processor.

Definition at line 987 of file dof_map.C.

References _send_list, Elem::active(), MeshBase::active_elements_begin(), MeshBase::active_elements_end(), Elem::active_family_tree_by_neighbor(), MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), dof_indices(), end_dof(), first_dof(), MeshBase::max_node_id(), Elem::n_neighbors(), Elem::n_nodes(), Elem::neighbor(), Elem::node(), DofObject::processor_id(), and libMesh::processor_id().

Referenced by distribute_dofs().

{
  START_LOG("add_neighbors_to_send_list()", "DofMap");
  
  //-------------------------------------------------------------------------
  // 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 LIBMESH_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 remote DOF indices into the send list
                  for (unsigned int j=0; j != di.size(); ++j)
                    if (di[j] < this->first_dof() ||
                        di[j] >= this->end_dof())
                      _send_list.push_back(di[j]);
                }
          }
    }

  // 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 remote DOF indices into the send list
            for (unsigned int j=0; j != di.size(); ++j)
              if (di[j] < this->first_dof() ||
                  di[j] >= this->end_dof())
                _send_list.push_back(di[j]);
        }
    }
  
  STOP_LOG("add_neighbors_to_send_list()", "DofMap");
}

void DofMap::add_periodic_boundary

(

const PeriodicBoundary &

periodic_boundary

)

Adds a copy of the specified periodic boundary to the system.

Definition at line 1304 of file dof_map_constraints.C.

References _periodic_boundaries, PeriodicBoundary::myboundary, PeriodicBoundary::pairedboundary, and PeriodicBoundary::translation_vector.

{
  PeriodicBoundary boundary = periodic_boundary;
  PeriodicBoundary inverse_boundary;
  inverse_boundary.myboundary = boundary.pairedboundary;
  inverse_boundary.pairedboundary = boundary.myboundary;
  inverse_boundary.translation_vector = -boundary.translation_vector;

  std::pair<unsigned int, PeriodicBoundary> bp
    (boundary.myboundary, boundary);
  std::pair<unsigned int, PeriodicBoundary> ibp
    (boundary.pairedboundary, inverse_boundary);

  _periodic_boundaries.insert(bp);
  _periodic_boundaries.insert(ibp);
}

void DofMap::add_variable

(

const System::Variable &

var

)

Add an unknown of order order and finite element type type to the system of equations.

Definition at line 56 of file dof_map.C.

References _variables.

{
  _variables.push_back (var);
}

void DofMap::allgather_recursive_constraints

(

)

Gathers any relevant constraint equations from other processors

Definition at line 881 of file dof_map_constraints.C.

References _dof_constraints, _end_df, is_constrained_dof(), std::max(), libMesh::n_processors(), and libMesh::processor_id().

Referenced by create_dof_constraints().

{
  // This function must be run on all processors at once
  parallel_only();

  // Return immediately if there's nothing to gather
  if (libMesh::n_processors() == 1)
    return;

  // We might get to return immediately if none of the processors
  // found any constraints
  unsigned int has_constraints = !_dof_constraints.empty();
  Parallel::max(has_constraints);
  if (!has_constraints)
    return;

  // We might have calculated constraints for constrained dofs
  // which live on other processors.
  // Push these out first.
  {
  std::vector<std::vector<unsigned int> > pushed_ids(libMesh::n_processors());
  std::vector<unsigned int> pushed_on_proc(libMesh::n_processors(), 0);

  // Count the constraints to push to each processor
  unsigned int push_proc_id = 0;
  for (DofConstraints::iterator i = _dof_constraints.begin();
         i != _dof_constraints.end(); ++i)
    {
      unsigned int constrained = i->first;
      while (constrained >= _end_df[push_proc_id])
        push_proc_id++;
      pushed_on_proc[push_proc_id]++;
    }
  for (unsigned int p = 0; p != libMesh::n_processors(); ++p)
    pushed_ids[p].reserve(pushed_on_proc[p]);

  // Collect the constraints to push to each processor
  push_proc_id = 0;
  for (DofConstraints::iterator i = _dof_constraints.begin();
         i != _dof_constraints.end(); ++i)
    {
      unsigned int constrained = i->first;
      while (constrained >= _end_df[push_proc_id])
        push_proc_id++;
      pushed_ids[push_proc_id].push_back(constrained);
    }

  // Now trade constraint rows
  for (unsigned int p = 0; p != libMesh::n_processors(); ++p)
    {
      // Push to processor procup while receiving from procdown
      unsigned int procup = (libMesh::processor_id() + p) %
                             libMesh::n_processors();
      unsigned int procdown = (libMesh::n_processors() +
                               libMesh::processor_id() - p) %
                               libMesh::n_processors();

      // Pack the constraint rows to push to procup
      std::vector<std::vector<unsigned int> > pushed_keys(pushed_ids[procup].size());
      std::vector<std::vector<Real> > pushed_vals(pushed_ids[procup].size());
      for (unsigned int i = 0; i != pushed_ids[procup].size(); ++i) 
        {
          DofConstraintRow &row = _dof_constraints[pushed_ids[procup][i]];
          unsigned int row_size = row.size();
          pushed_keys[i].reserve(row_size);
          pushed_vals[i].reserve(row_size);
          for (DofConstraintRow::iterator j = row.begin();
               j != row.end(); ++j)
            {
              pushed_keys[i].push_back(j->first);
              pushed_vals[i].push_back(j->second);
            }
        }

      // Trade pushed constraint rows
      std::vector<unsigned int> pushed_ids_to_me;
      std::vector<std::vector<unsigned int> > pushed_keys_to_me;
      std::vector<std::vector<Real> > pushed_vals_to_me;
      Parallel::send_receive(procup, pushed_ids[procup],
                             procdown, pushed_ids_to_me);
      Parallel::send_receive(procup, pushed_keys,
                             procdown, pushed_keys_to_me);
      Parallel::send_receive(procup, pushed_vals,
                             procdown, pushed_vals_to_me);
      libmesh_assert (pushed_ids_to_me.size() == pushed_keys_to_me.size());
      libmesh_assert (pushed_ids_to_me.size() == pushed_vals_to_me.size());

      // Add the constraints that I've been sent
      for (unsigned int i = 0; i != pushed_ids_to_me.size(); ++i)
        {
          libmesh_assert (pushed_keys_to_me[i].size() == pushed_vals_to_me[i].size());

          unsigned int constrained = pushed_ids_to_me[i];

          // If we don't already have a constraint for this dof,
          // add the one we were sent
          if (!this->is_constrained_dof(constrained))
            {
              DofConstraintRow &row = _dof_constraints[constrained];
              for (unsigned int j = 0; j != pushed_keys_to_me[i].size(); ++j)
                {
                  row[pushed_keys_to_me[i][j]] = pushed_vals_to_me[i][j];
                }
            }
        }
    }
  }

  // Now start checking for any other constraints we need
  // to know about, requesting them recursively.

  // Create a set containing the DOFs we already depend on
  typedef std::set<unsigned int> RCSet;
  RCSet unexpanded_set;

  for (DofConstraints::iterator i = _dof_constraints.begin();
         i != _dof_constraints.end(); ++i)
    {
      unexpanded_set.insert(i->first);
    }

  // We have to keep recursing while the unexpanded set is
  // nonempty on *any* processor
  unsigned int unexpanded_set_nonempty = !unexpanded_set.empty();
  Parallel::max(unexpanded_set_nonempty);

  while (unexpanded_set_nonempty)
    {
      // Request set
      RCSet request_set;

      // Request sets to send to each processor
      std::vector<std::vector<unsigned int> >
        requested_ids(libMesh::n_processors());

      // And the sizes of each
      std::vector<unsigned int> ids_on_proc(libMesh::n_processors(), 0);

      // Fill (and thereby sort and uniq!) the main request set
      for (RCSet::iterator i = unexpanded_set.begin();
           i != unexpanded_set.end(); ++i)
        {
          DofConstraintRow &row = _dof_constraints[*i];
          for (DofConstraintRow::iterator j = row.begin();
               j != row.end(); ++j)
            request_set.insert(j->first);
        }

      // Clear the unexpanded constraint set; we're about to expand it
      unexpanded_set.clear();

      // Count requests by processor
      unsigned int proc_id = 0;
      for (RCSet::iterator i = request_set.begin();
           i != request_set.end(); ++i)
        {
          while (*i >= _end_df[proc_id])
            proc_id++;
          ids_on_proc[proc_id]++;
        }
      for (unsigned int p = 0; p != libMesh::n_processors(); ++p)
        requested_ids[p].reserve(ids_on_proc[p]);

      // Prepare each processor's request set
      proc_id = 0;
      for (RCSet::iterator i = request_set.begin();
           i != request_set.end(); ++i)
        {
          while (*i >= _end_df[proc_id])
            proc_id++;
          requested_ids[proc_id].push_back(*i);
        }

      // Now request constraint rows from other processors
      for (unsigned int p=1; p != libMesh::n_processors(); ++p)
        {
          // Trade my requests with processor procup and procdown
          unsigned int procup = (libMesh::processor_id() + p) %
                                 libMesh::n_processors();
          unsigned int procdown = (libMesh::n_processors() +
                                   libMesh::processor_id() - p) %
                                   libMesh::n_processors();
          std::vector<unsigned int> request_to_fill;
          Parallel::send_receive(procup, requested_ids[procup],
                                 procdown, request_to_fill);

          // Fill those requests
          std::vector<std::vector<unsigned int> > row_keys(request_to_fill.size());
          std::vector<std::vector<Real> > row_vals(request_to_fill.size());
          for (unsigned int i=0; i != request_to_fill.size(); ++i)
            {
              unsigned int constrained = request_to_fill[i];
              if (_dof_constraints.count(constrained))
                {
                  DofConstraintRow &row = _dof_constraints[constrained];
                  unsigned int row_size = row.size();
                  row_keys[i].reserve(row_size);
                  row_vals[i].reserve(row_size);
                  for (DofConstraintRow::iterator j = row.begin();
                       j != row.end(); ++j)
                    {
                      row_keys[i].push_back(j->first);
                      row_vals[i].push_back(j->second);
                    }
                }
            }

          // Trade back the results
          std::vector<std::vector<unsigned int> > filled_keys;
          std::vector<std::vector<Real> > filled_vals;
          Parallel::send_receive(procdown, row_keys,
                                 procup, filled_keys);
          Parallel::send_receive(procdown, row_vals,
                                 procup, filled_vals);
          libmesh_assert (filled_keys.size() == requested_ids[procup].size());
          libmesh_assert (filled_vals.size() == requested_ids[procup].size());

          // Add any new constraint rows we've found
          for (unsigned int i=0; i != requested_ids[procup].size(); ++i)
            {
              libmesh_assert (filled_keys[i].size() == filled_vals[i].size());
              if (!filled_keys[i].empty())
                {
                  unsigned int constrained = requested_ids[procup][i];
                  DofConstraintRow &row = _dof_constraints[constrained];
                  for (unsigned int j = 0; j != filled_keys[i].size(); ++j)
                    row[filled_keys[i][j]] = filled_vals[i][j];

                  // And prepare to check for more recursive constraints
                  unexpanded_set.insert(constrained);
                }
            }
        }

      // We have to keep recursing while the unexpanded set is
      // nonempty on *any* processor
      unexpanded_set_nonempty = !unexpanded_set.empty();
      Parallel::max(unexpanded_set_nonempty);
    }
}

void DofMap::attach_matrix

(

SparseMatrix< Number > &

matrix

)

Additional matrices may be handled with this DofMap. They are initialized to the same sparsity structure as the major matrix.

Definition at line 90 of file dof_map.C.

References _matrices, and SparseMatrix< T >::attach_dof_map().

Referenced by EigenSystem::init_data(), and ImplicitSystem::init_matrices().

{
  _matrices.push_back(&matrix);
  
  matrix.attach_dof_map (*this);
}

void DofMap::build_constraint_matrix	(	DenseMatrix< Number > &	C,
		std::vector< unsigned int > &	elem_dofs,
		const bool	called_recursively = false
	)			const [private]

Build the constraint matrix C associated with the element degree of freedom indices elem_dofs. The optional parameter called_recursively should be left at the default value false. This is used to handle the special case of an element's degrees of freedom being constrained in terms of other, local degrees of freedom. The usual case is for an elements DOFs to be constrained by some other, external DOFs.

Definition at line 761 of file dof_map_constraints.C.

References _dof_constraints, is_constrained_dof(), DenseMatrixBase< T >::m(), DenseMatrixBase< T >::n(), DenseMatrix< T >::resize(), and DenseMatrix< T >::right_multiply().

Referenced by constrain_element_dyad_matrix(), constrain_element_matrix(), constrain_element_matrix_and_vector(), constrain_element_vector(), constrain_nothing(), extract_local_vector(), and max_constraint_error().

{
  if (!called_recursively) START_LOG("build_constraint_matrix()", "DofMap");

  // Create a set containing the DOFs we already depend on
  typedef std::set<unsigned int> RCSet;
  RCSet dof_set;

  bool we_have_constraints = false;

  // Next insert any other dofs the current dofs 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]))
      {
        we_have_constraints = true;

        // 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;
        
// Constraint rows in p refinement may be empty
//      libmesh_assert (!constraint_row.empty());
        
        for (DofConstraintRow::const_iterator
               it=constraint_row.begin(); it != constraint_row.end();
             ++it)
          dof_set.insert (it->first);
      }

  // May be safe to return at this point
  // (but remember to stop the perflog)
  if (!we_have_constraints)
    {
      STOP_LOG("build_constraint_matrix()", "DofMap");
      return;
    }

  // delay inserting elem_dofs for efficiency in the case of
  // no constraints.  In that case we don't get here!
  dof_set.insert (elem_dofs.begin(),
                  elem_dofs.end());

  // If we added any DOFS then we need to do this recursively.
  // It is possible that we just added a DOF that is also
  // constrained!
  //
  // Also, we need to handle the special case of an element having DOFs
  // constrained in terms of other, local DOFs
  if ((dof_set.size() != elem_dofs.size()) || // case 1: constrained in terms of other DOFs
      !called_recursively)                    // case 2: constrained in terms of our own DOFs
    {
      // Create a new list of element DOFs containing the
      // contents of the current dof_set.
      std::vector<unsigned int> new_elem_dofs (dof_set.begin(),
                                               dof_set.end());
      
      // Now we can build the constraint matrix.
      // Note that resize also zeros for a DenseMatrix<Number>.
      C.resize (elem_dofs.size(), new_elem_dofs.size());
      
      // Create the C constraint matrix.
      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;
            
// p refinement creates empty constraint rows
//          libmesh_assert (!constraint_row.empty());
            
            for (DofConstraintRow::const_iterator
                   it=constraint_row.begin(); it != constraint_row.end();
                 ++it)
              for (unsigned int j=0; j<new_elem_dofs.size(); j++)
                if (new_elem_dofs[j] == it->first)
                  C(i,j) = it->second;
          }
        else
          {
            for (unsigned int j=0; j<new_elem_dofs.size(); j++)
              if (new_elem_dofs[j] == elem_dofs[i])
                C(i,j) =  1.;
          }     

      // May need to do this recursively.  It is possible
      // that we just replaced a constrained DOF with another
      // constrained DOF.
      elem_dofs = new_elem_dofs;
      
      DenseMatrix<Number> Cnew;
      
      this->build_constraint_matrix (Cnew, elem_dofs, true);

      if ((C.n() == Cnew.m()) &&
          (Cnew.n() == elem_dofs.size())) // If the constraint matrix                                    
        C.right_multiply(Cnew);           // is constrained...
      
      libmesh_assert (C.n() == elem_dofs.size());
    }
  
  if (!called_recursively) STOP_LOG("build_constraint_matrix()", "DofMap");  
}

void DofMap::clear

(

)

Free all memory associated with the object, but keep the mesh pointer.

Definition at line 622 of file dof_map.C.

References _dof_constraints, _end_df, _first_df, _matrices, _n_dfs, _n_nz, _n_old_dfs, _n_oz, _send_list, _var_first_local_df, and _variables.

Referenced by ~DofMap().

{
  // we don't want to clear
  // the coupling matrix!
  // It should not change...
  //_dof_coupling->clear();

  _variables.clear();  
  _first_df.clear();
  _end_df.clear();
  _var_first_local_df.clear();
  _send_list.clear();
  _n_nz.clear();
  _n_oz.clear();


#ifdef LIBMESH_ENABLE_AMR

  _dof_constraints.clear();
  _n_old_dfs = 0;

#endif

  _matrices.clear();

  _n_dfs = 0;
}

void DofMap::compute_sparsity

(

const MeshBase &

mesh

)

Computes the sparsity pattern for the matrix corresponding to proc_id. Produces data that can be fed to Petsc for preallocation of sparse matrices.

Definition at line 1138 of file dof_map.C.

References _dof_coupling, _matrices, _n_nz, _n_oz, MeshBase::active_elements_begin(), MeshBase::active_elements_end(), MeshBase::is_prepared(), n_dofs_on_processor(), SparsityPattern::Build::n_nz, SparsityPattern::Build::n_oz, n_variables(), Threads::parallel_reduce(), MeshBase::processor_id(), SparsityPattern::Build::sparsity_pattern, and use_coupled_neighbor_dofs().

Referenced by EigenSystem::init_data(), and ImplicitSystem::init_matrices().

{
  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)

  // 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;

  // See if we need to include sparsity pattern entries for coupling
  // between neighbor dofs
  bool implicit_neighbor_dofs = this->use_coupled_neighbor_dofs(mesh);
  
  // 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::constrain_element_dyad_matrix	(	DenseVector< Number > &	v,
		DenseVector< Number > &	w,
		std::vector< unsigned int > &	row_dofs,
		bool	asymmetric_constraint_rows = true
	)			const [inline]

Constrains a dyadic element matrix B = v w'. This method requires the element matrix to be square, in which case the elem_dofs correspond to the global DOF indices of both the rows and columns of the element matrix. For this case the rows and columns of the matrix necessarily correspond to variables of the same approximation order.

Definition at line 520 of file dof_map_constraints.C.

References build_constraint_matrix(), is_constrained_dof(), DenseMatrixBase< T >::m(), DenseMatrixBase< T >::n(), n_constrained_dofs(), DenseVector< T >::resize(), DenseVector< T >::size(), and DenseMatrix< T >::transpose().

{
  libmesh_assert (v.size() == row_dofs.size());
  libmesh_assert (w.size() == row_dofs.size());

  // check for easy return
  if (this->n_constrained_dofs() == 0)
    return;
  
  // The constrained RHS is built up as R^T F.  
  DenseMatrix<Number> R;

  this->build_constraint_matrix (R, row_dofs);

  START_LOG("cnstrn_elem_dyad_mat()", "DofMap");
    
  // It is possible that the vector is not constrained at all.
  if ((R.m() == v.size()) &&
      (R.n() == row_dofs.size())) // if the RHS is constrained
    {
      // Compute the matrix-vector products
      DenseVector<Number> old_v(v);
      DenseVector<Number> old_w(w);

      // resize RHS & zero before summation
      v.resize(row_dofs.size());
      w.resize(row_dofs.size());

      // compute matrix/vector product
      for (unsigned int i=0; i<row_dofs.size(); i++)
        for (unsigned int j=0; j<old_v.size(); j++)
          {
            v(i) += R.transpose(i,j)*old_v(j);
            w(i) += R.transpose(i,j)*old_w(j);
          }
      
      libmesh_assert (row_dofs.size() == v.size());
      libmesh_assert (row_dofs.size() == w.size());

      /* Constrain only v, not w.  */

      for (unsigned int i=0; i<row_dofs.size(); i++)
        if (this->is_constrained_dof(row_dofs[i]))
          {     
            // If the DOF is constrained
            v(i) = 0.;
          }
    } // end if the RHS is constrained.
  
  STOP_LOG("cnstrn_elem_dyad_mat()", "DofMap");
}

void DofMap::constrain_element_matrix	(	DenseMatrix< Number > &	matrix,
		std::vector< unsigned int > &	row_dofs,
		std::vector< unsigned int > &	col_dofs,
		bool	asymmetric_constraint_rows = true
	)			const [inline]

Constrains the element matrix. This method allows the element matrix to be non-square, in which case the row_dofs and col_dofs may be of different size and correspond to variables approximated in different spaces.

Definition at line 391 of file dof_map_constraints.C.

References _dof_constraints, build_constraint_matrix(), is_constrained_dof(), DenseMatrix< T >::left_multiply_transpose(), DenseMatrixBase< T >::m(), DenseMatrixBase< T >::n(), n_constrained_dofs(), and DenseMatrix< T >::right_multiply().

{
  libmesh_assert (row_dofs.size() == matrix.m());
  libmesh_assert (col_dofs.size() == matrix.n());

  // check for easy return
  if (this->n_constrained_dofs() == 0)
    return;
  
  // The constrained matrix is built up as R^T K C.
  DenseMatrix<Number> R;
  DenseMatrix<Number> C;

  // Safeguard against the user passing us the same
  // object for row_dofs and col_dofs.  If that is done
  // the calls to build_matrix would fail
  std::vector<unsigned int> orig_row_dofs(row_dofs);
  std::vector<unsigned int> orig_col_dofs(col_dofs);
  
  this->build_constraint_matrix (R, orig_row_dofs);
  this->build_constraint_matrix (C, orig_col_dofs);

  START_LOG("constrain_elem_matrix()", "DofMap");
  
  row_dofs = orig_row_dofs;
  col_dofs = orig_col_dofs;
  
    
  // It is possible that the matrix is not constrained at all.
  if ((R.m() == matrix.m()) &&
      (R.n() == row_dofs.size()) &&
      (C.m() == matrix.n()) &&
      (C.n() == col_dofs.size())) // If the matrix is constrained
    {
      // K_constrained = R^T K C
      matrix.left_multiply_transpose  (R);
      matrix.right_multiply (C);
      
      
      libmesh_assert (matrix.m() == row_dofs.size());
      libmesh_assert (matrix.n() == col_dofs.size());
      
      
      for (unsigned int i=0; i<row_dofs.size(); i++)
        if (this->is_constrained_dof(row_dofs[i]))
          {
            for (unsigned int j=0; j<matrix.n(); j++)
              matrix(i,j) = 0.;
          
            // If the DOF is constrained
            matrix(i,i) = 1.;
            
            if (asymmetric_constraint_rows)
              {
                DofConstraints::const_iterator
                  pos = _dof_constraints.find(row_dofs[i]);
            
                libmesh_assert (pos != _dof_constraints.end());
            
                const DofConstraintRow& constraint_row = pos->second;
            
                libmesh_assert (!constraint_row.empty());
            
                for (DofConstraintRow::const_iterator
                       it=constraint_row.begin(); it != constraint_row.end();
                     ++it)
                  for (unsigned int j=0; j<col_dofs.size(); j++)
                    if (col_dofs[j] == it->first)
                      matrix(i,j) = -it->second;        
              }
          }
    } // end if is constrained...
  
  STOP_LOG("constrain_elem_matrix()", "DofMap");  
}

void DofMap::constrain_element_matrix	(	DenseMatrix< Number > &	matrix,
		std::vector< unsigned int > &	elem_dofs,
		bool	asymmetric_constraint_rows = true
	)			const [inline]

Constrains the element matrix. This method requires the element matrix to be square, in which case the elem_dofs correspond to the global DOF indices of both the rows and columns of the element matrix. For this case the rows and columns of the matrix necessarily correspond to variables of the same approximation order.

If asymmetric_constraint_rows is set to true (as it is by default), constraint row equations will be reinforced in a way which breaks matrix symmetry but makes inexact linear solver solutions more likely to satisfy hanging node constraints.

Definition at line 226 of file dof_map_constraints.C.

References _dof_constraints, build_constraint_matrix(), is_constrained_dof(), DenseMatrix< T >::left_multiply_transpose(), DenseMatrixBase< T >::m(), DenseMatrixBase< T >::n(), n_constrained_dofs(), and DenseMatrix< T >::right_multiply().

Referenced by FEMSystem::assembly().

{
  libmesh_assert (elem_dofs.size() == matrix.m());
  libmesh_assert (elem_dofs.size() == matrix.n());

  // check for easy return
  if (this->n_constrained_dofs() == 0)
    return;
    
  // The constrained matrix is built up as C^T K C.    
  DenseMatrix<Number> C;

  
  this->build_constraint_matrix (C, elem_dofs);

  START_LOG("constrain_elem_matrix()", "DofMap");
  
  // It is possible that the matrix is not constrained at all.
  if ((C.m() == matrix.m()) &&
      (C.n() == elem_dofs.size())) // It the matrix is constrained
    {
      // Compute the matrix-matrix-matrix product C^T K C
      matrix.left_multiply_transpose  (C);
      matrix.right_multiply (C);
      
      
      libmesh_assert (matrix.m() == matrix.n());
      libmesh_assert (matrix.m() == elem_dofs.size());
      libmesh_assert (matrix.n() == elem_dofs.size());
      
      
      for (unsigned int i=0; i<elem_dofs.size(); i++)
        // If the DOF is constrained
        if (this->is_constrained_dof(elem_dofs[i]))
          {
            for (unsigned int j=0; j<matrix.n(); j++)
              matrix(i,j) = 0.;
            
            matrix(i,i) = 1.;
            
            if (asymmetric_constraint_rows)
              {
                DofConstraints::const_iterator
                  pos = _dof_constraints.find(elem_dofs[i]);
            
                libmesh_assert (pos != _dof_constraints.end());
            
                const DofConstraintRow& constraint_row = pos->second;
            
                libmesh_assert (!constraint_row.empty());
            
                for (DofConstraintRow::const_iterator
                       it=constraint_row.begin(); it != constraint_row.end();
                     ++it)
                  for (unsigned int j=0; j<elem_dofs.size(); j++)
                    if (elem_dofs[j] == it->first)
                      matrix(i,j) = -it->second;        
              }
          }
    } // end if is constrained...
  
  STOP_LOG("constrain_elem_matrix()", "DofMap");  
}

void DofMap::constrain_element_matrix_and_vector	(	DenseMatrix< Number > &	matrix,
		DenseVector< Number > &	rhs,
		std::vector< unsigned int > &	elem_dofs,
		bool	asymmetric_constraint_rows = true
	)			const [inline]

Constrains the element matrix and vector. This method requires the element matrix to be square, in which case the elem_dofs correspond to the global DOF indices of both the rows and columns of the element matrix. For this case the rows and columns of the matrix necessarily correspond to variables of the same approximation order.

Definition at line 294 of file dof_map_constraints.C.

References _dof_constraints, build_constraint_matrix(), is_constrained_dof(), DenseMatrix< T >::left_multiply_transpose(), DenseMatrixBase< T >::m(), DenseMatrixBase< T >::n(), n_constrained_dofs(), DenseVector< T >::resize(), DenseMatrix< T >::right_multiply(), DenseVector< T >::size(), DenseMatrix< T >::transpose(), and DenseVector< T >::zero().

Referenced by FEMSystem::assembly().

{
  libmesh_assert (elem_dofs.size() == matrix.m());
  libmesh_assert (elem_dofs.size() == matrix.n());
  libmesh_assert (elem_dofs.size() == rhs.size());

  // check for easy return
  if (this->n_constrained_dofs() == 0)
    return;
  
  // The constrained matrix is built up as C^T K C.
  // The constrained RHS is built up as C^T F
  DenseMatrix<Number> C;
  
  this->build_constraint_matrix (C, elem_dofs);
  
  START_LOG("cnstrn_elem_mat_vec()", "DofMap");
  
  // It is possible that the matrix is not constrained at all.
  if ((C.m() == matrix.m()) &&
      (C.n() == elem_dofs.size())) // It the matrix is constrained
    {
      // Compute the matrix-matrix-matrix product C^T K C
      matrix.left_multiply_transpose  (C);
      matrix.right_multiply (C);
      
      
      libmesh_assert (matrix.m() == matrix.n());
      libmesh_assert (matrix.m() == elem_dofs.size());
      libmesh_assert (matrix.n() == elem_dofs.size());
      

      for (unsigned int i=0; i<elem_dofs.size(); i++)
        if (this->is_constrained_dof(elem_dofs[i]))
          {
            for (unsigned int j=0; j<matrix.n(); j++)
              matrix(i,j) = 0.;
            
            // If the DOF is constrained
            matrix(i,i) = 1.;
            
            // This will put a nonsymmetric entry in the constraint
            // row to ensure that the linear system produces the
            // correct value for the constrained DOF.
            if (asymmetric_constraint_rows)
              {
                DofConstraints::const_iterator
                  pos = _dof_constraints.find(elem_dofs[i]);
            
                libmesh_assert (pos != _dof_constraints.end());
            
                const DofConstraintRow& constraint_row = pos->second;
            
// p refinement creates empty constraint rows
//          libmesh_assert (!constraint_row.empty());
            
                for (DofConstraintRow::const_iterator
                       it=constraint_row.begin(); it != constraint_row.end();
                     ++it)
                  for (unsigned int j=0; j<elem_dofs.size(); j++)
                    if (elem_dofs[j] == it->first)
                      matrix(i,j) = -it->second;        
              }
          }

      
      // Compute the matrix-vector product C^T F
      DenseVector<Number> old_rhs(rhs);

      // resize the RHS vector & 0 before summation
      rhs.resize(elem_dofs.size());
      rhs.zero();
      
      // compute matrix/vector product
      for (unsigned int i=0; i<elem_dofs.size(); i++)
        for (unsigned int j=0; j<old_rhs.size(); j++)
          rhs(i) += C.transpose(i,j)*old_rhs(j);
        

      libmesh_assert (elem_dofs.size() == rhs.size());

      for (unsigned int i=0; i<elem_dofs.size(); i++)
        if (this->is_constrained_dof(elem_dofs[i]))
          {     
            // If the DOF is constrained
            rhs(i) = 0.;
          }
    } // end if is constrained...
  
  STOP_LOG("cnstrn_elem_mat_vec()", "DofMap");  
}

void DofMap::constrain_element_vector	(	DenseVector< Number > &	rhs,
		std::vector< unsigned int > &	dofs,
		bool	asymmetric_constraint_rows = true
	)			const [inline]

Constrains the element vector.

Definition at line 472 of file dof_map_constraints.C.

References build_constraint_matrix(), is_constrained_dof(), DenseMatrixBase< T >::m(), DenseMatrixBase< T >::n(), n_constrained_dofs(), DenseVector< T >::resize(), DenseVector< T >::size(), DenseMatrix< T >::transpose(), and DenseVector< T >::zero().

Referenced by FEMSystem::assemble_qoi_derivative(), and FEMSystem::assembly().

{
  libmesh_assert (rhs.size() == row_dofs.size());

  // check for easy return
  if (this->n_constrained_dofs() == 0)
    return;
  
  // The constrained RHS is built up as R^T F.  
  DenseMatrix<Number> R;

  this->build_constraint_matrix (R, row_dofs);

  START_LOG("constrain_elem_vector()", "DofMap");
    
  // It is possible that the vector is not constrained at all.
  if ((R.m() == rhs.size()) &&
      (R.n() == row_dofs.size())) // if the RHS is constrained
    {
      // Compute the matrix-vector product
      DenseVector<Number> old_rhs(rhs);

      // resize RHS & zero before summation
      rhs.resize(row_dofs.size());
      rhs.zero();

      // compute matrix/vector product
      for (unsigned int i=0; i<row_dofs.size(); i++)
        for (unsigned int j=0; j<old_rhs.size(); j++)
          rhs(i) += R.transpose(i,j)*old_rhs(j);
      
      libmesh_assert (row_dofs.size() == rhs.size());

      for (unsigned int i=0; i<row_dofs.size(); i++)
        if (this->is_constrained_dof(row_dofs[i]))
          {     
            // If the DOF is constrained
            rhs(i) = 0.;
          }
    } // end if the RHS is constrained.
  
  STOP_LOG("constrain_elem_vector()", "DofMap");  
}

void DofMap::constrain_nothing

(

std::vector< unsigned int > &

dofs

)

const

Does not actually constrain anything, but modifies dofs in the same way as any of the constrain functions would do, i.e. adds those dofs in terms of which any of the existing dofs is constrained.

Definition at line 577 of file dof_map_constraints.C.

References build_constraint_matrix(), and n_constrained_dofs().

{
  // check for easy return
  if (this->n_constrained_dofs() == 0)
    return;

  // All the work is done by \p build_constraint_matrix.  We just need
  // a dummy matrix.
  DenseMatrix<Number> R;
  this->build_constraint_matrix (R, dofs);
}

void DofMap::constrain_p_dofs	(	unsigned int	var,
		const Elem *	elem,
		unsigned int	s,
		unsigned int	p
	)

Constrains degrees of freedom on side s of element elem which correspond to variable number var and to p refinement levels above p.

Definition at line 1242 of file dof_map_constraints.C.

References _dof_constraints, Elem::dim(), dim, DofObject::dof_number(), Elem::get_node(), Elem::is_node_on_side(), Elem::is_vertex(), DofObject::n_comp(), FEInterface::n_dofs_at_node(), Elem::n_nodes(), n_nodes, FEType::order, Elem::p_level(), Threads::spin_mtx, sys_number(), Elem::type(), type, and variable_type().

Referenced by FEBase::compute_periodic_constraints(), and FEBase::compute_proj_constraints().

{
  // We're constraining dofs on elem which correspond to p refinement
  // levels above p - this only makes sense if elem's p refinement
  // level is above p.
  libmesh_assert(elem->p_level() > p);
  libmesh_assert(s < elem->n_sides());

  const unsigned int sys_num = this->sys_number();
  const unsigned int dim = elem->dim();
  ElemType type = elem->type();
  FEType low_p_fe_type = this->variable_type(var);
  FEType high_p_fe_type = this->variable_type(var);
  low_p_fe_type.order = static_cast<Order>(low_p_fe_type.order + p);
  high_p_fe_type.order = static_cast<Order>(high_p_fe_type.order + 
                                            elem->p_level());

  const unsigned int n_nodes = elem->n_nodes();
  for (unsigned int n = 0; n != n_nodes; ++n)
    if (elem->is_node_on_side(n, s))
      {
        const Node * const node = elem->get_node(n);
        const unsigned int low_nc =
          FEInterface::n_dofs_at_node (dim, low_p_fe_type, type, n);
        const unsigned int high_nc =
          FEInterface::n_dofs_at_node (dim, high_p_fe_type, type, n);
        
        // since we may be running this method concurretly 
        // on multiple threads we need to acquire a lock 
        // before modifying the _dof_constraints object.
        Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
        
        if (elem->is_vertex(n))
          {
            // Add "this is zero" constraint rows for high p vertex
            // dofs
            for (unsigned int i = low_nc; i != high_nc; ++i)
              _dof_constraints[node->dof_number(sys_num,var,i)].clear();
          }
        else
          {
            const unsigned int total_dofs = node->n_comp(sys_num, var);
            libmesh_assert(total_dofs >= high_nc);
            // Add "this is zero" constraint rows for high p
            // non-vertex dofs, which are numbered in reverse
            for (unsigned int j = low_nc; j != high_nc; ++j)
              {
                const unsigned int i = total_dofs - j - 1;
                _dof_constraints[node->dof_number(sys_num,var,i)].clear();
              }
          }
      }
}

void DofMap::create_dof_constraints

(

const MeshBase &

mesh

)

Rebuilds the raw degree of freedom constraints.

Definition at line 112 of file dof_map_constraints.C.

References _dof_constraints, _periodic_boundaries, allgather_recursive_constraints(), dim, MeshBase::elements_begin(), MeshBase::elements_end(), MeshBase::is_prepared(), MeshBase::is_serial(), MeshBase::local_elements_begin(), MeshBase::local_elements_end(), MeshBase::mesh_dimension(), n_variables(), Threads::parallel_for(), and StoredRange< iterator_type, object_type >::reset().

Referenced by EquationSystems::reinit().

{
  START_LOG("create_dof_constraints()", "DofMap");

  libmesh_assert (mesh.is_prepared());
  
  const unsigned int dim = mesh.mesh_dimension();

  // Constraints are not necessary in 1D
  if (dim == 1)
  {
    // make sure we stop logging though
    STOP_LOG("create_dof_constraints()", "DofMap");
    return;
  }
  
  // Here we build the hanging node constraints.  This is done
  // by enforcing the condition u_a = u_b along hanging sides.
  // u_a = u_b is collocated at the nodes of side a, which gives
  // one row of the constraint matrix.
  
  // clear any existing constraints.
  _dof_constraints.clear();

  // define the range of elements of interest
  ConstElemRange range;
  {
    // With SerialMesh or a serial ParallelMesh, every processor
    // computes every constraint
    MeshBase::const_element_iterator
      elem_begin = mesh.elements_begin(),
      elem_end   = mesh.elements_end();
    
    // With a parallel ParallelMesh, processors compute only
    // their local constraints
    if (!mesh.is_serial())
      {
        elem_begin = mesh.local_elements_begin();
        elem_end   = mesh.local_elements_end(); 
      }

    // set the range to contain the specified elements
    range.reset (elem_begin, elem_end);
  }
  
  // Look at all the variables in the system.  Reset the element
  // range at each iteration -- there is no need to reconstruct it.
  for (unsigned int variable_number=0; variable_number<this->n_variables();
       ++variable_number, range.reset())
    Threads::parallel_for (range,
                           ComputeConstraints (_dof_constraints,
                                               *this,
#ifdef LIBMESH_ENABLE_PERIODIC
                                               _periodic_boundaries,
#endif
                                               mesh,
                                               variable_number));

  // With a parallelized Mesh, we've computed our local constraints,
  // but they may depend on non-local constraints that we'll need to
  // take into account.
  if (!mesh.is_serial())
    this->allgather_recursive_constraints();
  
  STOP_LOG("create_dof_constraints()", "DofMap");
}

void DofMap::distribute_dofs

(

MeshBase &

mesh

)

Distrubute dofs on the current mesh. Also builds the send list for processor proc_id, which defaults to 0 for ease of use in serial applications.

Definition at line 652 of file dof_map.C.

References _end_df, _first_df, _n_dfs, _n_old_dfs, _send_list, add_neighbors_to_send_list(), distribute_local_dofs_node_major(), distribute_local_dofs_var_major(), elem_ptr(), MeshBase::elements_begin(), MeshBase::elements_end(), invalidate_dofs(), MeshBase::is_prepared(), libMesh::n_processors(), node_ptr(), MeshBase::nodes_begin(), MeshBase::nodes_end(), libMesh::on_command_line(), libMesh::processor_id(), reinit(), and set_nonlocal_dof_objects().

Referenced by EquationSystems::reinit().

{
  // This function must be run on all processors at once
  parallel_only();

  // Log how long it takes to distribute the degrees of freedom
  START_LOG("distribute_dofs()", "DofMap");

  libmesh_assert (mesh.is_prepared());

  const unsigned int proc_id = libMesh::processor_id();
  const unsigned int n_proc  = libMesh::n_processors();
  
//  libmesh_assert (this->n_variables() > 0);
  libmesh_assert (proc_id < n_proc);
  
  // re-init in case the mesh has changed
  this->reinit(mesh);
  
  // By default distribute variables in a
  // var-major fashion, but allow run-time
  // specification
  bool node_major_dofs = libMesh::on_command_line ("--node_major_dofs");

  // The DOF counter, will be incremented as we encounter
  // new degrees of freedom
  unsigned int next_free_dof = 0;

  // Clear the send list before we rebuild it
  _send_list.clear();

  // Set temporary DOF indices on this processor
  if (node_major_dofs)
    this->distribute_local_dofs_node_major (next_free_dof, mesh);
  else
    this->distribute_local_dofs_var_major (next_free_dof, mesh);

  // Get DOF counts on all processors
  std::vector<unsigned int> dofs_on_proc(n_proc, 0);
  Parallel::allgather(next_free_dof, dofs_on_proc);

  // Resize the _first_df and _end_df arrays
  _first_df.resize(n_proc);
  _end_df.resize (n_proc);

  // Get DOF offsets
  _first_df[0] = 0;
  for (unsigned int i=1; i < n_proc; ++i)
    _first_df[i] = _end_df[i-1] = _first_df[i-1] + dofs_on_proc[i-1];
  _end_df[n_proc-1] = _first_df[n_proc-1] + dofs_on_proc[n_proc-1];

  // Clear all the current DOF indices
  // (distribute_dofs expects them cleared!)
  this->invalidate_dofs(mesh);

  next_free_dof = _first_df[proc_id];

  // Set permanent DOF indices on this processor
  if (node_major_dofs)
    this->distribute_local_dofs_node_major (next_free_dof, mesh);
  else
    this->distribute_local_dofs_var_major (next_free_dof, mesh);

  libmesh_assert(next_free_dof == _end_df[proc_id]);

  //------------------------------------------------------------
  // At this point, all n_comp and dof_number values on local
  // DofObjects should be correct, but a ParallelMesh might have
  // incorrect values on non-local DofObjects.  Let's request the
  // correct values from each other processor.

  if (libMesh::n_processors() > 1)
    {
      this->set_nonlocal_dof_objects(mesh.nodes_begin(),
                                     mesh.nodes_end(),
                                     mesh, &DofMap::node_ptr);

      this->set_nonlocal_dof_objects(mesh.elements_begin(),
                                     mesh.elements_end(),
                                     mesh, &DofMap::elem_ptr);
    }
  
  // Set the total number of degrees of freedom
#ifdef LIBMESH_ENABLE_AMR
  _n_old_dfs = _n_dfs;
#endif
  _n_dfs = _end_df[n_proc-1];

  STOP_LOG("distribute_dofs()", "DofMap");

  // Note that in the add_neighbors_to_send_list nodes on processor
  // boundaries that are shared by multiple elements are added for
  // each element.
  this->add_neighbors_to_send_list(mesh);
  
  // Here we used to clean up that data structure; now System and
  // EquationSystems call that for us, after we've added constraint
  // dependencies to the send_list too.
  // this->sort_send_list ();
}

void DofMap::distribute_local_dofs_node_major	(	unsigned int &	next_free_dof,
		MeshBase &	mesh
	)			[private]

Distributes the global degrees of freedom, for dofs on this processor. In this format all the degrees of freedom at a node/element are in contiguous blocks. Note in particular that the degrees of freedom for a given variable are not in contiguous blocks, as in the case of distribute_local_dofs_var_major. Starts at index next_free_dof, and increments it to the post-final index. If build_send_list is true, builds the send list. If false, clears and reserves the send list

Definition at line 754 of file dof_map.C.

References _n_SCALAR_dofs, _var_first_local_df, MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), DofObject::dof_number(), first_dof(), Elem::get_node(), DofObject::invalid_id, MeshTools::libmesh_assert_valid_node_procids(), MeshBase::local_nodes_begin(), MeshBase::local_nodes_end(), DofObject::n_comp(), Elem::n_nodes(), n_nodes, libMesh::n_processors(), n_variables(), DofObject::n_vars(), n_vars, FEType::order, libMesh::processor_id(), DofObject::processor_id(), libMeshEnums::SCALAR, DofObject::set_dof_number(), Elem::subdomain_id(), sys_number(), System::Variable::type(), type, and variable().

Referenced by distribute_dofs().

{
  const unsigned int sys_num = this->sys_number();
  const unsigned int n_vars  = this->n_variables();

  // We now only add remote dofs to the _send_list
  // unsigned int send_list_size = 0;

  // _var_first_local_df does not work with node_major dofs
  _var_first_local_df.resize(n_vars+1);
  std::fill (_var_first_local_df.begin(),
             _var_first_local_df.end(),
             DofObject::invalid_id);

  //-------------------------------------------------------------------------
  // First count and assign temporary numbers to local dofs
  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);

          for (unsigned var=0; var<n_vars; var++)
          {
              if( (this->variable(var).type().family != SCALAR) &&
                  (this->variable(var).active_on_subdomain(elem->subdomain_id())) )
              {
                // 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);
                  }
              }
          }
        }

      // Now number the element DOFS
      for (unsigned var=0; var<n_vars; var++)
        if ( (this->variable(var).type().family != SCALAR) &&
             (this->variable(var).active_on_subdomain(elem->subdomain_id())) )
          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);
            }
    } // done looping over elements

  // Finally, count up the SCALAR dofs
  this->_n_SCALAR_dofs = 0;
  for (unsigned var=0; var<n_vars; var++)
    {
      if( this->variable(var).type().family == SCALAR )
        {
          this->_n_SCALAR_dofs += this->variable(var).type().order;
          continue;
        }
    }

  // Only increment next_free_dof if we're on the processor
  // that holds this SCALAR variable
  if ( libMesh::processor_id() == (libMesh::n_processors()-1) )
    next_free_dof += _n_SCALAR_dofs;

#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
}

void DofMap::distribute_local_dofs_var_major	(	unsigned int &	next_free_dof,
		MeshBase &	mesh
	)			[private]

Distributes the global degrees of freedom, for dofs on this processor. In this format the local degrees of freedom are in a contiguous block for each variable in the system. Starts at index next_free_dof, and increments it to the post-final index.

Definition at line 867 of file dof_map.C.

References _n_SCALAR_dofs, _var_first_local_df, MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), System::Variable::active_on_subdomain(), DofObject::dof_number(), FEType::family, first_dof(), Elem::get_node(), DofObject::invalid_id, MeshTools::libmesh_assert_valid_node_procids(), MeshBase::local_nodes_begin(), MeshBase::local_nodes_end(), DofObject::n_comp(), Elem::n_nodes(), n_nodes, libMesh::n_processors(), n_variables(), DofObject::n_vars(), n_vars, FEType::order, libMesh::processor_id(), DofObject::processor_id(), libMeshEnums::SCALAR, DofObject::set_dof_number(), Elem::subdomain_id(), sys_number(), type, System::Variable::type(), and variable().

Referenced by distribute_dofs().

{
  const unsigned int sys_num = this->sys_number();
  const unsigned int n_vars  = this->n_variables();

  // We now only add remote dofs to the _send_list
  // unsigned int send_list_size = 0;

  // We will cache the first local index for each variable
  _var_first_local_df.clear();

  //-------------------------------------------------------------------------
  // First count and assign temporary numbers to local dofs
  for (unsigned var=0; var<n_vars; var++)
    {
      _var_first_local_df.push_back(next_free_dof);

      const System::Variable var_description = this->variable(var);

      // Skip the SCALAR dofs
      if(var_description.type().family == SCALAR)
        continue;

      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;

          // ... and only variables which are active on
          // on this element's subdomain
          if (!var_description.active_on_subdomain(elem->subdomain_id()))
            continue;

          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);
                }
            }
                  
          // 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);
            }
        } // end loop on elements
    } // end loop on variables

  // Finally, count up the SCALAR dofs
  this->_n_SCALAR_dofs = 0;
  for (unsigned var=0; var<n_vars; var++)
    {
      if( this->variable(var).type().family == SCALAR )
        {
          this->_n_SCALAR_dofs += this->variable(var).type().order;
          continue;
        }
    }

  // Only increment next_free_dof if we're on the processor
  // that holds this SCALAR variable
  if ( libMesh::processor_id() == (libMesh::n_processors()-1) )
    next_free_dof += _n_SCALAR_dofs;

  // Cache the last local dof number too
  _var_first_local_df.push_back(next_free_dof);

#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
}

void DofMap::dof_indices	(	const Elem *const	elem,
		std::vector< unsigned int > &	di,
		const unsigned int	vn = libMesh::invalid_uint
	)			const

Fills the vector di with the global degree of freedom indices for the element. If no variable number is specified then all variables are returned.

Definition at line 1280 of file dof_map.C.

References Elem::active(), Elem::dim(), dim, DofObject::dof_number(), FEInterface::extra_hanging_dofs(), FEType::family, Elem::get_node(), DofObject::invalid_id, libMesh::invalid_uint, Elem::is_vertex(), libMeshEnums::LAGRANGE, DofObject::n_comp(), n_dofs(), FEInterface::n_dofs_at_node(), FEInterface::n_dofs_per_elem(), Elem::n_nodes(), n_nodes, DofObject::n_systems(), n_variables(), n_vars, FEType::order, Elem::p_level(), libMeshEnums::SCALAR, SCALAR_dof_indices(), Elem::subdomain_id(), sys_number(), Elem::type(), type, variable(), and variable_type().

Referenced by ExactSolution::_compute_error(), add_neighbors_to_send_list(), HPCoarsenTest::add_projection(), EquationSystems::build_discontinuous_solution_vector(), EquationSystems::build_solution_vector(), System::calculate_norm(), FEBase::coarsened_dof_values(), FEBase::compute_periodic_constraints(), FEBase::compute_proj_constraints(), JumpErrorEstimator::estimate_error(), ExactErrorEstimator::estimate_error(), MeshFunction::gradient(), MeshFunction::hessian(), System::local_dof_indices(), max_constraint_error(), FEMSystem::mesh_position_get(), System::ProjectVector::operator()(), PatchRecoveryErrorEstimator::EstimateError::operator()(), MeshFunction::operator()(), SparsityPattern::Build::operator()(), ErrorVector::plot_error(), FEMContext::reinit(), HPCoarsenTest::select_refinement(), EnsightIO::write_scalar_ascii(), and EnsightIO::write_vector_ascii().

{
  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

  // Create a vector to indicate which
  // SCALAR variables have been requested
  std::vector<unsigned int> SCALAR_var_numbers;
  SCALAR_var_numbers.clear();

  // Get the dof numbers
  for (unsigned int v=0; v<n_vars; v++)
    if ((v == vn) || (vn == libMesh::invalid_uint))
    {
      if(this->variable(v).type().family == SCALAR)
      {
        // We asked for this variable, so add it to the vector.
        SCALAR_var_numbers.push_back(v);

#ifdef DEBUG
        FEType fe_type = this->variable_type(v);
        tot_size += FEInterface::n_dofs(dim,
                                        fe_type,
                                        type);
#endif
      }
      else
      if (this->variable(v).active_on_subdomain(elem->subdomain_id()))
        { // 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);
                }
            }
        }
      } // end loop over variables

  // Finally append any SCALAR dofs that we asked for.
  std::vector<unsigned int> di_new;
  std::vector<unsigned int>::iterator it           = SCALAR_var_numbers.begin();
  std::vector<unsigned int>::const_iterator it_end = SCALAR_var_numbers.end();
  for( ; it != it_end; ++it)
  {
    this->SCALAR_dof_indices(di_new,*it);
    di.insert( di.end(), di_new.begin(), di_new.end());
  }

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

DofObject * DofMap::elem_ptr	(	MeshBase &	mesh,
		unsigned int	i
	)			const [private]

An adapter function that returns Elem pointers by index

Definition at line 106 of file dof_map.C.

References MeshBase::elem().

Referenced by distribute_dofs().

{
  return mesh.elem(i);
}

unsigned int DofMap::end_dof

(

const unsigned int

proc = libMesh::processor_id()

)

const [inline]

Returns the first dof index that is after all indices local to subdomain proc. Analogous to the end() member function of STL containers.

Definition at line 378 of file dof_map.h.

References _end_df.

Referenced by add_constraints_to_send_list(), add_neighbors_to_send_list(), enforce_constraints_exactly(), SparsityPattern::Build::join(), System::local_dof_indices(), SparsityPattern::Build::operator()(), and SparseMatrix< T >::print().

  { libmesh_assert(proc < _end_df.size()); return _end_df[proc]; }

void DofMap::enforce_constraints_exactly	(	const System &	system,
		NumericVector< Number > *	v = NULL
	)			const [inline]

Constrains the numeric vector v, which represents a solution defined on the mesh. This may need to be used after a linear solve, if your linear solver's solutions do not satisfy your DoF constraints to a tight enough tolerance.

If v == NULL, the system solution vector is constrained

Definition at line 591 of file dof_map_constraints.C.

References _dof_constraints, NumericVector< T >::close(), NumericVector< T >::closed(), end_dof(), first_dof(), AutoPtr< Tp >::get(), System::get_dof_map(), libMeshEnums::GHOSTED, NumericVector< T >::localize(), std::max(), n_constrained_dofs(), n_dofs(), n_local_dofs(), libMeshEnums::PARALLEL, libMeshEnums::SERIAL, NumericVector< T >::set(), NumericVector< T >::size(), System::solution, and NumericVector< T >::type().

Referenced by __libmesh_petsc_diff_solver_jacobian(), __libmesh_petsc_diff_solver_residual(), ImplicitSystem::adjoint_solve(), System::project_vector(), ImplicitSystem::sensitivity_solve(), PetscDiffSolver::solve(), NewtonSolver::solve(), ImplicitSystem::weighted_sensitivity_adjoint_solve(), and ImplicitSystem::weighted_sensitivity_solve().

{
  parallel_only();

  unsigned int local_constraints = this->n_constrained_dofs();
  Parallel::max(local_constraints);
  if (!local_constraints)
    return;
  
  START_LOG("enforce_constraints_exactly()","DofMap");

  if (!v)
    v = system.solution.get();

  NumericVector<Number> *v_local  = NULL; // will be initialized below
  NumericVector<Number> *v_global = NULL; // will be initialized below
  AutoPtr<NumericVector<Number> > v_built;
  if (v->type() == SERIAL)
    {
      v_built = NumericVector<Number>::build();
      v_built->init(this->n_dofs(), this->n_local_dofs(), true, PARALLEL);
      v_built->close();

      for (unsigned int i=v_built->first_local_index();
           i<v_built->last_local_index(); i++)
        v_built->set(i, (*v)(i));
      v_built->close();
      v_global = v_built.get();

      v_local = v;
      libmesh_assert (v_local->closed());
    }
  else if (v->type() == PARALLEL)
    {
      v_built = NumericVector<Number>::build();
      v_built->init (v->size(), v->size(), true, SERIAL);
      v->localize(*v_built);
      v_built->close();
      v_local = v_built.get();

      v_global = v;
    }
  else if (v->type() == GHOSTED)
    {
      v_local = v;
      v_global = v;
    }
  else // unknown v->type()
    {
      std::cerr << "ERROR: Unknown v->type() == " << v->type() 
                << std::endl;
      libmesh_error();
    }

  // We should never hit these asserts because we should error-out in
  // else clause above.  Just to be sure we don't try to use v_local
  // and v_global uninitialized...
  libmesh_assert(v_local);
  libmesh_assert(v_global);
  libmesh_assert (this == &(system.get_dof_map()));

  DofConstraints::const_iterator c_it = _dof_constraints.begin();
  const DofConstraints::const_iterator c_end = _dof_constraints.end();

  for ( ; c_it != c_end; ++c_it)
    {
      unsigned int constrained_dof = c_it->first;
      if (constrained_dof < this->first_dof() ||
          constrained_dof >= this->end_dof())
        continue;

      const DofConstraintRow constraint_row = c_it->second;

      Number exact_value = 0;
      for (DofConstraintRow::const_iterator
           j=constraint_row.begin(); j != constraint_row.end();
           ++j)
        exact_value += j->second * (*v_local)(j->first);

      v_global->set(constrained_dof, exact_value);
    }

  // If the old vector was serial, we probably need to send our values
  // to other processors
  if (v->type() == SERIAL)
    {
      v_global->localize (*v);
    }
  v->close();

  STOP_LOG("enforce_constraints_exactly()","DofMap");
}

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

Builds the local element vector Ue from the global vector Ug, accounting for any constrained degrees of freedom. For an element without constrained degrees of freedom this is the trivial mapping

Note that the user must ensure that the element vector Ue is properly sized when calling this method. This is because there is no resize() method in the DenseVectorBase<> class.

Definition at line 1210 of file dof_map.C.

References build_constraint_matrix(), DenseVectorBase< T >::el(), NumericVector< T >::first_local_index(), is_constrained_dof(), NumericVector< T >::last_local_index(), DenseMatrixBase< T >::m(), DenseMatrixBase< T >::n(), DenseVectorBase< T >::size(), and DenseVectorBase< T >::zero().

{
#ifdef LIBMESH_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()));

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

void DofMap::find_connected_dofs

(

std::vector< unsigned int > &

elem_dofs

)

const [private]

Finds all the DOFS associated with the element DOFs elem_dofs. This will account for off-element couplings via hanging nodes.

Definition at line 1737 of file dof_map.C.

References _dof_constraints, and is_constrained_dof().

Referenced by SparsityPattern::Build::operator()().

{
  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();      
      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)
}

unsigned int DofMap::first_dof

(

const unsigned int

proc = libMesh::processor_id()

)

const [inline]

Returns the first dof index that is local to subdomain proc.

Definition at line 363 of file dof_map.h.

References _first_df.

Referenced by add_constraints_to_send_list(), add_neighbors_to_send_list(), distribute_local_dofs_node_major(), distribute_local_dofs_var_major(), enforce_constraints_exactly(), SparsityPattern::Build::join(), System::local_dof_indices(), SparsityPattern::Build::operator()(), SparseMatrix< T >::print(), and set_nonlocal_dof_objects().

  { libmesh_assert(proc < _first_df.size()); return _first_df[proc]; }

std::string ReferenceCounter::get_info

(

)

[static, inherited]

Gets a string containing the reference information.

Definition at line 45 of file reference_counter.C.

References ReferenceCounter::_counts, and QuadratureRules::name().

Referenced by ReferenceCounter::print_info().

{
#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)

  std::ostringstream out;
  
  out << '\n'
      << " ---------------------------------------------------------------------------- \n"
      << "| Reference count information                                                |\n"
      << " ---------------------------------------------------------------------------- \n";
  
  for (Counts::iterator it = _counts.begin();
       it != _counts.end(); ++it)
    {
      const std::string name(it->first);
      const unsigned int creations    = it->second.first;
      const unsigned int destructions = it->second.second;

      out << "| " << name << " reference count information:\n"
          << "|  Creations:    " << creations    << '\n'
          << "|  Destructions: " << destructions << '\n';
    }
  
  out << " ---------------------------------------------------------------------------- \n";

  return out.str();

#else

  return "";
  
#endif
}

const std::vector<unsigned int>& DofMap::get_n_nz

(

)

const [inline]

Returns a constant reference to the _n_nz list for this processor. The vector contains the bandwidth of the on-processor coupling for each row of the global matrix that the current processor owns. This information can be used to preallocate space for a parallel sparse matrix.

Definition at line 299 of file dof_map.h.

References _n_nz.

Referenced by PetscMatrix< T >::init(), LaspackMatrix< T >::init(), and EpetraMatrix< T >::update_sparsity_pattern().

{ return _n_nz; }

const std::vector<unsigned int>& DofMap::get_n_oz

(

)

const [inline]

Returns a constant reference to the _n_oz list for this processor. The vector contains the bandwidth of the off-processor coupling for each row of the global matrix that the current processor owns. This information can be used to preallocate space for a parallel sparse matrix.

Definition at line 307 of file dof_map.h.

References _n_oz.

Referenced by PetscMatrix< T >::init(), LaspackMatrix< T >::init(), and EpetraMatrix< T >::update_sparsity_pattern().

{ return _n_oz; }

const std::vector<unsigned int>& DofMap::get_send_list

(

)

const [inline]

Returns a constant reference to the _send_list for this processor. The _send_list contains the global indices of all the variables in the global solution vector that influence the current processor. This information can be used for gathers at each solution step to retrieve solution values needed for computation.

Definition at line 291 of file dof_map.h.

References _send_list.

Referenced by UnsteadySolver::solve().

{ return _send_list; }

void ReferenceCounter::increment_constructor_count

(

const std::string &

name

)

[inline, protected, inherited]

Increments the construction counter. Should be called in the constructor of any derived class that will be reference counted.

Definition at line 149 of file reference_counter.h.

References ReferenceCounter::_counts, and Threads::spin_mtx.

Referenced by ReferenceCountedObject< SparseMatrix< T > >::ReferenceCountedObject().

{
  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
  std::pair<unsigned int, unsigned int>& p = _counts[name];

  p.first++;
}

void ReferenceCounter::increment_destructor_count

(

const std::string &

name

)

[inline, protected, inherited]

Increments the destruction counter. Should be called in the destructor of any derived class that will be reference counted.

Definition at line 167 of file reference_counter.h.

References ReferenceCounter::_counts, and Threads::spin_mtx.

Referenced by ReferenceCountedObject< SparseMatrix< T > >::~ReferenceCountedObject().

{
  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
  std::pair<unsigned int, unsigned int>& p = _counts[name];

  p.second++;
}

void DofMap::invalidate_dofs

(

MeshBase &

mesh

)

const [private]

Invalidates all active DofObject dofs for this system

Definition at line 601 of file dof_map.C.

References MeshBase::active_elements_begin(), MeshBase::active_elements_end(), MeshBase::nodes_begin(), MeshBase::nodes_end(), and sys_number().

Referenced by distribute_dofs(), and reinit().

{
  const unsigned int sys_num = this->sys_number();

  // 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)->invalidate_dofs(sys_num);
  
  // All the elements
  MeshBase::element_iterator       elem_it  = mesh.active_elements_begin();
  const MeshBase::element_iterator elem_end = mesh.active_elements_end(); 

  for ( ; elem_it != elem_end; ++elem_it)
    (*elem_it)->invalidate_dofs(sys_num);
}

bool DofMap::is_constrained_dof

(

const unsigned int

dof

)

const [inline]

Returns:

true if the degree of freedom dof is constrained, false otherwise.

Definition at line 895 of file dof_map.h.

References _dof_constraints.

Referenced by add_constraint_row(), allgather_recursive_constraints(), build_constraint_matrix(), FEBase::compute_periodic_constraints(), constrain_element_dyad_matrix(), constrain_element_matrix(), constrain_element_matrix_and_vector(), constrain_element_vector(), extract_local_vector(), find_connected_dofs(), max_constraint_error(), and process_constraints().

{
  if (_dof_constraints.count(dof) != 0)
    return true;

  return false;
}

bool DofMap::is_periodic_boundary

(

const unsigned int

boundaryid

)

const [inline]

Returns:

true if the boundary given by boundaryid is periodic, false otherwise

Definition at line 909 of file dof_map.h.

References _periodic_boundaries.

{
  if (_periodic_boundaries.count(boundaryid) != 0)
    return true;

  return false;
}

unsigned int DofMap::last_dof

(

const unsigned int

proc = libMesh::processor_id()

)

const [inline]

Returns the last dof index that is local to subdomain proc. This function is now deprecated, because it returns nonsense in the rare case where proc has no local dof indices. Use end_dof() instead.

Definition at line 371 of file dof_map.h.

References _end_df.

  { libmesh_deprecated(); libmesh_assert(proc < _end_df.size()); return (_end_df[proc] - 1); }  

std::pair< Real, Real > DofMap::max_constraint_error	(	const System &	system,
		NumericVector< Number > *	v = NULL
	)			const

Tests the constrained degrees of freedom on the numeric vector v, which represents a solution defined on the mesh, returning a pair whose first entry is the maximum absolute error on a constrained DoF and whose second entry is the maximum relative error. Useful for debugging purposes.

If v == NULL, the system solution vector is tested.

Definition at line 688 of file dof_map_constraints.C.

References MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), build_constraint_matrix(), NumericVector< T >::closed(), dof_indices(), NumericVector< T >::first_local_index(), System::get_dof_map(), System::get_mesh(), is_constrained_dof(), NumericVector< T >::last_local_index(), DenseMatrixBase< T >::m(), std::max(), mesh, DenseMatrixBase< T >::n(), and System::solution.

{
  if (!v)
    v = system.solution.get();
  NumericVector<Number> &vec = *v;

  // We'll assume the vector is closed
  libmesh_assert (vec.closed());

  Real max_absolute_error = 0., max_relative_error = 0.;

  const MeshBase &mesh = system.get_mesh();

  libmesh_assert (this == &(system.get_dof_map()));

  // indices on each element
  std::vector<unsigned int> local_dof_indices;

  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* elem = *elem_it;

      this->dof_indices(elem, local_dof_indices);
      std::vector<unsigned int> raw_dof_indices = local_dof_indices;

      // Constraint matrix for each element
      DenseMatrix<Number> C;

      this->build_constraint_matrix (C, local_dof_indices);

      // Continue if the element is unconstrained
      if (!C.m())
        continue;

      libmesh_assert(C.m() == raw_dof_indices.size());
      libmesh_assert(C.n() == local_dof_indices.size());

      for (unsigned int i=0; i!=C.m(); ++i)
        {
          // Recalculate any constrained dof owned by this processor
          unsigned int global_dof = raw_dof_indices[i];
          if (this->is_constrained_dof(global_dof) &&
              global_dof >= vec.first_local_index() &&
              global_dof < vec.last_local_index())
          {
            Number exact_value = 0;
            for (unsigned int j=0; j!=C.n(); ++j)
              {
                if (local_dof_indices[j] != global_dof)
                  exact_value += C(i,j) * 
                    vec(local_dof_indices[j]);
              }

            max_absolute_error = std::max(max_absolute_error,
              std::abs(vec(global_dof) - exact_value));
            max_relative_error = std::max(max_relative_error,
              std::abs(vec(global_dof) - exact_value)
              / std::abs(exact_value));
          }
        }
    }

  return std::pair<Real, Real>(max_absolute_error, max_relative_error);
}

unsigned int DofMap::n_constrained_dofs

(

)

const [inline]

Returns:

the total number of constrained degrees of freedom in the problem.

Definition at line 452 of file dof_map.h.

References _dof_constraints.

Referenced by constrain_element_dyad_matrix(), constrain_element_matrix(), constrain_element_matrix_and_vector(), constrain_element_vector(), constrain_nothing(), and enforce_constraints_exactly().

{ return _dof_constraints.size(); }

unsigned int DofMap::n_dofs

(

)

const [inline]

Returns:

the total number of degrees of freedom in the problem.

Definition at line 342 of file dof_map.h.

References _n_dfs.

Referenced by dof_indices(), enforce_constraints_exactly(), PetscMatrix< T >::init(), LaspackMatrix< T >::init(), SparsityPattern::Build::join(), old_dof_indices(), UnsteadySolver::old_nonlinear_solution(), SCALAR_dof_indices(), and EpetraMatrix< T >::update_sparsity_pattern().

{ return _n_dfs; }

unsigned int DofMap::n_dofs_on_processor

(

const unsigned int

proc

)

const [inline]

Returns the number of degrees of freedom on subdomain proc.

Definition at line 357 of file dof_map.h.

References _end_df, and _first_df.

Referenced by compute_sparsity(), PetscMatrix< T >::init(), LaspackMatrix< T >::init(), SparsityPattern::Build::join(), n_local_dofs(), SparsityPattern::Build::operator()(), and EpetraMatrix< T >::update_sparsity_pattern().

  { libmesh_assert(proc < _first_df.size()); return (_end_df[proc] - _first_df[proc]); }

unsigned int DofMap::n_local_dofs

(

)

const [inline]

Returns:

the number of degrees of freedom on this processor.

Definition at line 351 of file dof_map.h.

References n_dofs_on_processor(), and libMesh::processor_id().

Referenced by enforce_constraints_exactly().

  { return this->n_dofs_on_processor (libMesh::processor_id()); }

static unsigned int ReferenceCounter::n_objects

(

)

[inline, static, inherited]

Prints the number of outstanding (created, but not yet destroyed) objects.

Definition at line 76 of file reference_counter.h.

References ReferenceCounter::_n_objects.

  { return _n_objects; }

unsigned int DofMap::n_old_dofs

(

)

const [inline]

Returns:

the total number of degrees of freedom on old_dof_objects

Definition at line 640 of file dof_map.h.

References _n_old_dfs.

Referenced by SCALAR_dof_indices().

{ return _n_old_dfs; }

unsigned int DofMap::n_SCALAR_dofs

(

)

const [inline]

Returns:

the number of SCALAR dofs.

Definition at line 346 of file dof_map.h.

References _n_SCALAR_dofs.

Referenced by SCALAR_dof_indices().

{ return _n_SCALAR_dofs; }

unsigned int DofMap::n_variables

(

)

const [inline]

Returns the number of variables in the global solution vector. Defaults to 1, should be 1 for a scalar equation, 3 for 2D incompressible Navier Stokes (u,v,p), etc...

Definition at line 335 of file dof_map.h.

References _variables.

Referenced by compute_sparsity(), create_dof_constraints(), distribute_local_dofs_node_major(), distribute_local_dofs_var_major(), dof_indices(), old_dof_indices(), SparsityPattern::Build::operator()(), reinit(), SCALAR_dof_indices(), set_nonlocal_dof_objects(), and use_coupled_neighbor_dofs().

  { return _variables.size(); }

DofObject * DofMap::node_ptr	(	MeshBase &	mesh,
		unsigned int	i
	)			const [private]

An adapter function that returns Node pointers by index

Definition at line 99 of file dof_map.C.

References MeshBase::node_ptr().

Referenced by distribute_dofs().

{
  return mesh.node_ptr(i);
}

void DofMap::old_dof_indices	(	const Elem *const	elem,
		std::vector< unsigned int > &	di,
		const unsigned int	vn = libMesh::invalid_uint
	)			const

After a mesh is refined and repartitioned it is possible that the _send_list will need to be augmented. This is the case when an element is refined and its children end up on different processors than the parent. These children will need values from the parent when projecting the solution onto the refined mesh, hence the parent's DOF indices need to be included in the _send_list. Fills the vector di with the global degree of freedom indices for the element using the DofMap::old_dof_object. If no variable number is specified then all variables are returned.

Definition at line 1495 of file dof_map.C.

References Elem::active(), Elem::dim(), dim, DofObject::dof_number(), FEInterface::extra_hanging_dofs(), FEType::family, Elem::get_node(), DofObject::invalid_id, libMesh::invalid_uint, Elem::is_vertex(), Elem::JUST_COARSENED, Elem::JUST_REFINED, libMeshEnums::LAGRANGE, DofObject::n_comp(), n_dofs(), FEInterface::n_dofs_at_node(), FEInterface::n_dofs_per_elem(), Elem::n_nodes(), n_nodes, DofObject::n_systems(), n_variables(), n_vars, DofObject::old_dof_object, FEType::order, Elem::p_level(), Elem::p_refinement_flag(), Elem::refinement_flag(), libMeshEnums::SCALAR, SCALAR_dof_indices(), Elem::subdomain_id(), sys_number(), Elem::type(), type, variable(), and variable_type().

Referenced by FEBase::coarsened_dof_values(), and System::ProjectVector::operator()().

{
  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

  // Create a vector to indicate which
  // SCALAR variables have been requested
  std::vector<unsigned int> SCALAR_var_numbers;
  SCALAR_var_numbers.clear();

  // Get the dof numbers
  for (unsigned int v=0; v<n_vars; v++)
    if ((v == vn) || (vn == libMesh::invalid_uint))
    {
      if(this->variable(v).type().family == SCALAR)
        {
          // We asked for this variable, so add it to the vector.
          SCALAR_var_numbers.push_back(v);

#ifdef DEBUG
          FEType fe_type = this->variable_type(v);
          tot_size += FEInterface::n_dofs(dim,
                                          fe_type,
                                          type);
#endif
        }
      else
      if (this->variable(v).active_on_subdomain(elem->subdomain_id()))
        { // 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);
                }
            }
        }
    } // end loop over variables

  // Finally append any SCALAR dofs that we asked for.
  std::vector<unsigned int> di_new;
  std::vector<unsigned int>::iterator it           = SCALAR_var_numbers.begin();
  std::vector<unsigned int>::const_iterator it_end = SCALAR_var_numbers.end();
  for( ; it != it_end; ++it)
  {
    this->SCALAR_dof_indices(di_new,*it,true);
    di.insert( di.end(), di_new.begin(), di_new.end());
  }
  
#ifdef DEBUG
  libmesh_assert (tot_size == di.size());
#endif
  
  STOP_LOG("old_dof_indices()", "DofMap");  
}

void DofMap::prepare_send_list

(

)

Takes the _send_list vector (which may have duplicate entries) and sorts it. The duplicate entries are then removed, resulting in a sorted _send_list with unique entries.

Definition at line 1089 of file dof_map.C.

References _send_list.

Referenced by EquationSystems::reinit().

{
  START_LOG("prepare_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("prepare_send_list()", "DofMap");
}

void DofMap::print_dof_constraints

(

std::ostream &

os = std::cout

)

const

Prints the _dof_constraints data structure.

Definition at line 201 of file dof_map_constraints.C.

References _dof_constraints.

{
  os << "DOF CONSTRAINTS OUTPUT:"
     << std::endl;
  
  for (DofConstraints::const_iterator it=_dof_constraints.begin();
       it != _dof_constraints.end(); ++it)
    {
      const unsigned int i = it->first;
      const DofConstraintRow& row = it->second;

      os << "Constraints for DOF " << i
         << ": \t";

      for (DofConstraintRow::const_iterator pos=row.begin();
           pos != row.end(); ++pos)
        os << " (" << pos->first << ","
           << pos->second << ")\t";

      os << std::endl;
    }
}

void ReferenceCounter::print_info

(

)

[static, inherited]

Prints the reference information to std::cout.

Definition at line 83 of file reference_counter.C.

References ReferenceCounter::get_info().

{
#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
  
  std::cout << ReferenceCounter::get_info();
  
#endif
}

void DofMap::process_constraints

(

)

Postprocesses any constrained degrees of freedom in elem_dofs to be constrained only in terms of unconstrained dofs, then adds unconstrained dofs to the send_list and prepares that for use. This should be run after both system (create_dof_constraints) and user constraints have all been added.

Definition at line 1124 of file dof_map_constraints.C.

References _dof_constraints, add_constraints_to_send_list(), and is_constrained_dof().

Referenced by EquationSystems::reinit().

{
  // Create a set containing the DOFs we already depend on
  typedef std::set<unsigned int> RCSet;
  RCSet unexpanded_set;

  for (DofConstraints::iterator i = _dof_constraints.begin();
         i != _dof_constraints.end(); ++i)
    unexpanded_set.insert(i->first);

  while (!unexpanded_set.empty())
    for (RCSet::iterator i = unexpanded_set.begin();
         i != unexpanded_set.end(); /* nothing */)
      {
        // If the DOF is constrained
        DofConstraints::iterator
          pos = _dof_constraints.find(*i);
        
        libmesh_assert (pos != _dof_constraints.end());
        
        DofConstraintRow& constraint_row = pos->second;

        std::vector<unsigned int> constraints_to_expand;

        for (DofConstraintRow::const_iterator
               it=constraint_row.begin(); it != constraint_row.end();
             ++it)
          if (this->is_constrained_dof(it->first))
            {
              unexpanded_set.insert(it->first);
              constraints_to_expand.push_back(it->first);
            }

        for (unsigned int j = 0; j != constraints_to_expand.size();
             ++j)
          {
            unsigned int expandable = constraints_to_expand[j];

            DofConstraints::const_iterator
              subpos = _dof_constraints.find(expandable);
        
            libmesh_assert (subpos != _dof_constraints.end());
        
            const DofConstraintRow& subconstraint_row = subpos->second;
            
            for (DofConstraintRow::const_iterator
                   it=subconstraint_row.begin();
                   it != subconstraint_row.end(); ++it)
              {
                constraint_row[it->first] += it->second *
                                constraint_row[expandable];
              }
            constraint_row.erase(expandable);
          }

        if (constraints_to_expand.empty())
          unexpanded_set.erase(i++);
        else
          i++;
      }

  // Now that we have our root constraint dependencies sorted out, add
  // them to the send_list
  this->add_constraints_to_send_list();
}

void DofMap::reinit

(

MeshBase &

mesh

)

Reinitialize the underlying data strucures conformal to the current mesh.

Definition at line 263 of file dof_map.C.

References _n_SCALAR_dofs, MeshBase::active_elements_begin(), MeshBase::active_elements_end(), System::Variable::active_on_subdomain(), Elem::dim(), dim, DofObject::dof_number(), MeshBase::elements_begin(), MeshBase::elements_end(), Utility::enum_to_string(), FEInterface::extra_hanging_dofs(), FEType::family, Elem::get_node(), DofObject::has_dofs(), invalidate_dofs(), MeshBase::is_prepared(), Elem::is_vertex(), Elem::JUST_REFINED, std::max(), FEInterface::max_order(), DofObject::n_comp(), FEInterface::n_dofs_at_node(), FEInterface::n_dofs_per_elem(), Elem::n_nodes(), n_variables(), MeshBase::nodes_begin(), MeshBase::nodes_end(), DofObject::old_dof_object, FEType::order, Elem::p_level(), Elem::refinement_flag(), libMeshEnums::SCALAR, DofObject::set_dof_number(), DofObject::set_n_comp(), DofObject::set_old_dof_object(), Elem::set_p_level(), Elem::subdomain_id(), sys_number(), Elem::type(), type, variable(), and variable_type().

Referenced by distribute_dofs().

{
  libmesh_assert (mesh.is_prepared());
  
  START_LOG("reinit()", "DofMap");
  
  //this->clear();

  const unsigned int n_var = this->n_variables();

#ifdef LIBMESH_ENABLE_AMR
  
  //------------------------------------------------------------
  // Clear the old_dof_objects for all the nodes
  // and elements so that we can overwrite them
  {
    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)->clear_old_dof_object();
        libmesh_assert ((*node_it)->old_dof_object == NULL);
      }
    
    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)->clear_old_dof_object();
        libmesh_assert ((*elem_it)->old_dof_object == NULL);
      }
  }

  
  //------------------------------------------------------------
  // Set the old_dof_objects for the elements that
  // weren't just created, if these old dof objects
  // had variables
  {
    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* elem = *elem_it;

        // Skip the elements that were just refined
        if (elem->refinement_flag() == Elem::JUST_REFINED) continue;

        for (unsigned int n=0; n<elem->n_nodes(); n++)
          {
            Node* node = elem->get_node(n);

            if (node->old_dof_object == NULL)
              if (node->has_dofs())
                node->set_old_dof_object();
          }

        libmesh_assert (elem->old_dof_object == NULL);

        if (elem->has_dofs())
          elem->set_old_dof_object();
      }
  }

#endif // #ifdef LIBMESH_ENABLE_AMR

  
  //------------------------------------------------------------
  // Then set the number of variables for each \p DofObject
  // equal to n_variables() for this system.  This will
  // handle new \p DofObjects that may have just been created
  {
    // 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)->set_n_vars(this->sys_number(),n_var);

    // 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)->set_n_vars(this->sys_number(),n_var);
  }


  // Zero _n_SCALAR_dofs, it will be updated below.
  this->_n_SCALAR_dofs = 0;

  //------------------------------------------------------------
  // Next allocate space for the DOF indices
  for (unsigned int var=0; var<this->n_variables(); var++)
    {
      const System::Variable &var_description = this->variable(var);
      const FEType& base_fe_type              = this->variable_type(var);

      // Don't need to loop over elements for a SCALAR variable
      // Just increment _n_SCALAR_dofs
      if(base_fe_type.family == SCALAR)
      {
        this->_n_SCALAR_dofs += base_fe_type.order;
        continue;
      }

      // This should be constant even on p-refined elements
      const bool extra_hanging_dofs =
        FEInterface::extra_hanging_dofs(base_fe_type);

      // For all the active elements
      MeshBase::element_iterator       elem_it  = mesh.active_elements_begin();
      const MeshBase::element_iterator elem_end = mesh.active_elements_end(); 
 
      // Count vertex degrees of freedom first
      for ( ; elem_it != elem_end; ++elem_it)
        {
          Elem* elem  = *elem_it;
          libmesh_assert (elem != NULL);
          
          // Skip the numbering if this variable is 
          // not active on this element's subdoman
          if (!var_description.active_on_subdomain(elem->subdomain_id()))
            continue;

          const ElemType type = elem->type();
          const unsigned int dim = elem->dim();

          FEType fe_type = base_fe_type;

#ifdef LIBMESH_ENABLE_AMR
          // Make sure we haven't done more p refinement than we can
          // handle
          if (elem->p_level() + base_fe_type.order >
              FEInterface::max_order(base_fe_type, type))
            {
#  ifdef DEBUG
              if (FEInterface::max_order(base_fe_type,type) <
                  static_cast<unsigned int>(base_fe_type.order))
                {
                  std::cerr << "ERROR: Finite element "
                    << Utility::enum_to_string(base_fe_type.family)
                    << " on geometric element "
                    << Utility::enum_to_string(type) << std::endl
                    << "only supports FEInterface::max_order = " 
                    << FEInterface::max_order(base_fe_type,type)
                    << ", not fe_type.order = " << base_fe_type.order
                    << std::endl;

                  libmesh_error();
                }

              std::cerr << "WARNING: Finite element "
                    << Utility::enum_to_string(base_fe_type.family)
                    << " on geometric element "
                    << Utility::enum_to_string(type) << std::endl
                    << "could not be p refined past FEInterface::max_order = " 
                    << FEInterface::max_order(base_fe_type,type)
                    << std::endl;
#  endif
              elem->set_p_level(FEInterface::max_order(base_fe_type,type)
                                - base_fe_type.order);
            }
#endif

          fe_type.order = static_cast<Order>(fe_type.order +
                                             elem->p_level());
             
          // Allocate the vertex DOFs
          for (unsigned int n=0; n<elem->n_nodes(); n++)
            {
              Node* node = elem->get_node(n);

              if (elem->is_vertex(n))
                {
                  const unsigned int old_node_dofs =
                    node->n_comp(this->sys_number(), var);

                  const unsigned int vertex_dofs =
                    std::max(FEInterface::n_dofs_at_node(dim, fe_type,
                                                         type, n),
                             old_node_dofs);
                  
                  // Some discontinuous FEs have no vertex dofs
                  if (vertex_dofs > old_node_dofs)
                    {
                      node->set_n_comp(this->sys_number(), var,
                                       vertex_dofs);
                      // Abusing dof_number to set a "this is a
                      // vertex" flag
                      node->set_dof_number(this->sys_number(),
                                           var, 0, vertex_dofs);
                    }
                }
            }
        } // done counting vertex dofs

      // count edge & face dofs next
      elem_it = mesh.active_elements_begin();

      for ( ; elem_it != elem_end; ++elem_it)
        {
          Elem* elem = *elem_it;
          libmesh_assert (elem != NULL);
          
          // Skip the numbering if this variable is 
          // not active on this element's subdoman
          if (!var_description.active_on_subdomain(elem->subdomain_id()))
            continue;

          const ElemType type = elem->type();
          const unsigned int dim = elem->dim();
             
          FEType fe_type = base_fe_type;
          fe_type.order = static_cast<Order>(fe_type.order +
                                             elem->p_level());

          // Allocate the edge and face DOFs
          for (unsigned int n=0; n<elem->n_nodes(); n++)
            {
              Node* node = elem->get_node(n);

              const unsigned int old_node_dofs =
                node->n_comp(this->sys_number(), var);

              const unsigned int vertex_dofs = old_node_dofs?
                node->dof_number (this->sys_number(),var,0):0;

              const unsigned int new_node_dofs =
                FEInterface::n_dofs_at_node(dim, fe_type, type, n);

              // We've already allocated vertex DOFs
              if (elem->is_vertex(n))
                {
                  libmesh_assert(old_node_dofs >= vertex_dofs);
                  libmesh_assert(vertex_dofs >= new_node_dofs);
                }
              // We need to allocate the rest
              else
                {
                  // If this has no dofs yet, it needs no vertex
                  // dofs, so we just give it edge or face dofs
                  if (!old_node_dofs)
                    {
                      node->set_n_comp(this->sys_number(), var,
                                       new_node_dofs);
                      // Abusing dof_number to set a "this has no
                      // vertex dofs" flag
                      if (new_node_dofs)
                        node->set_dof_number(this->sys_number(),
                                             var, 0, 0);
                    }

                  // If this has dofs, but has no vertex dofs,
                  // it may still need more edge or face dofs if
                  // we're p-refined.
                  else if (vertex_dofs == 0)
                    {
                      if (new_node_dofs > old_node_dofs)
                        {
                          node->set_n_comp(this->sys_number(), var,
                                           new_node_dofs);
                          node->set_dof_number(this->sys_number(),
                                               var, 0, vertex_dofs);
                        }
                    }
                  // If this is another element's vertex,
                  // add more (non-overlapping) edge/face dofs if
                  // necessary
                  else if (extra_hanging_dofs)
                    {
                      if (new_node_dofs > old_node_dofs - vertex_dofs)
                        {
                          node->set_n_comp(this->sys_number(), var,
                                           vertex_dofs + new_node_dofs);
                          node->set_dof_number(this->sys_number(),
                                               var, 0, vertex_dofs);
                        }
                    }
                  // If this is another element's vertex, add any
                  // (overlapping) edge/face dofs if necessary
                  else
                    {
                      libmesh_assert(old_node_dofs >= vertex_dofs);
                      if (new_node_dofs > old_node_dofs)
                        {
                          node->set_n_comp(this->sys_number(), var,
                                           new_node_dofs);
                          node->set_dof_number(this->sys_number(),
                                               var, 0, vertex_dofs);
                        }
                    }
                }
            }
          // Allocate the element DOFs
          const unsigned int dofs_per_elem =
                          FEInterface::n_dofs_per_elem(dim, fe_type,
                                                       type);

          elem->set_n_comp(this->sys_number(), var, dofs_per_elem);

        }
    }

  // Calling DofMap::reinit() by itself makes little sense,
  // so we won't bother with nonlocal DofObjects.
  // Those will be fixed by distribute_dofs
/*
  //------------------------------------------------------------
  // At this point, all n_comp and dof_number values on local
  // DofObjects should be correct, but a ParallelMesh might have
  // incorrect values on non-local DofObjects.  Let's request the
  // correct values from each other processor.

  this->set_nonlocal_n_comps(mesh.nodes_begin(),
                             mesh.nodes_end(),
                             mesh,
                             &DofMap::node_ptr);

  this->set_nonlocal_n_comps(mesh.elements_begin(),
                             mesh.elements_end(),
                             mesh,
                             &DofMap::elem_ptr);
*/

  //------------------------------------------------------------
  // Finally, clear all the current DOF indices
  // (distribute_dofs expects them cleared!)
  this->invalidate_dofs(mesh);
  
  STOP_LOG("reinit()", "DofMap");
}

void DofMap::SCALAR_dof_indices	(	std::vector< unsigned int > &	di,
		const unsigned int	vn,
		const bool	old_dofs = false
	)			const

Fills the vector di with the global degree of freedom indices corresponding to the SCALAR variable vn. If old_dofs=true, the old SCALAR dof indices are returned. Note that we do not need to pass in an element since SCALARs are global variables.

Definition at line 1438 of file dof_map.C.

References FEType::family, n_dofs(), n_old_dofs(), n_SCALAR_dofs(), n_variables(), FEType::order, libMeshEnums::SCALAR, System::Variable::type(), type, and variable().

Referenced by dof_indices(), old_dof_indices(), and System::project_vector().

{
  START_LOG("SCALAR_dof_indices()", "DofMap");

  if(this->variable(vn).type().family != SCALAR)
  {
    std::cerr << "ERROR: SCALAR_dof_indices called for a non-SCALAR variable."
              << std::endl;
  }

  // Clear the DOF indices vector
  di.clear();

  // First we need to find out the first dof
  // index for each SCALAR.
  unsigned int first_SCALAR_dof_index = (old_dofs ? n_old_dofs() : n_dofs()) - n_SCALAR_dofs();
  std::map<unsigned int, unsigned int> SCALAR_first_dof_index;
  SCALAR_first_dof_index.clear();

  // Iterate over _all_ of the SCALARs and store each one's first dof index
  // We need to do this since the SCALAR dofs are packed contiguously
  for (unsigned int v=0; v<this->n_variables(); v++)
    if(this->variable(v).type().family == SCALAR)
    {
      unsigned int current_n_SCALAR_dofs = this->variable(v).type().order;
      SCALAR_first_dof_index.insert(
        std::pair<unsigned int, unsigned int>(v,first_SCALAR_dof_index) );
      first_SCALAR_dof_index += current_n_SCALAR_dofs;
    }

  // Now use vn to index into SCALAR_first_dof_index
  std::map<unsigned int, unsigned int>::const_iterator iter =
    SCALAR_first_dof_index.find(vn);

#ifdef DEBUG
  libmesh_assert(iter != SCALAR_first_dof_index.end());
#endif

  unsigned int current_first_SCALAR_dof_index = iter->second;

  // Also, get the number of SCALAR dofs from the variable order
  unsigned int current_n_SCALAR_dofs = this->variable(vn).type().order;

  for(unsigned int j=0; j<current_n_SCALAR_dofs; j++)
  {
    unsigned int index = current_first_SCALAR_dof_index+j;
    di.push_back(index);
  }

  STOP_LOG("SCALAR_dof_indices()", "DofMap");
}

template<typename iterator_type >

void DofMap::set_nonlocal_dof_objects	(	iterator_type	objects_begin,
		iterator_type	objects_end,
		MeshBase &	mesh,
		dofobject_accessor	objects
	)			[inline, private]

Helper function for distributing dofs in parallel

Definition at line 114 of file dof_map.C.

References _send_list, DofObject::dof_number(), first_dof(), DofObject::id(), DofObject::invalid_id, DofObject::invalid_processor_id, DofObject::n_comp(), libMesh::n_processors(), n_variables(), DofObject::n_vars(), libMesh::processor_id(), DofObject::processor_id(), DofObject::set_dof_number(), DofObject::set_n_comp(), and sys_number().

Referenced by distribute_dofs().

{
  // This function must be run on all processors at once
  parallel_only();

  // First, iterate over local objects to find out how many
  // are on each processor
  std::vector<unsigned int>
    ghost_objects_from_proc(libMesh::n_processors(), 0);

  iterator_type it  = objects_begin;

  for (; it != objects_end; ++it)
    {
      DofObject *obj = *it;

      if (obj)
        {
          unsigned int obj_procid = obj->processor_id();
          // We'd better be completely partitioned by now
          libmesh_assert(obj_procid != DofObject::invalid_processor_id);
          ghost_objects_from_proc[obj_procid]++;
        }
    }

  std::vector<unsigned int> objects_on_proc(libMesh::n_processors(), 0);
  Parallel::allgather(ghost_objects_from_proc[libMesh::processor_id()],
                      objects_on_proc);

#ifdef DEBUG
  for (unsigned int p=0; p != libMesh::n_processors(); ++p)
    libmesh_assert(ghost_objects_from_proc[p] <= objects_on_proc[p]);
#endif

  // Request sets to send to each processor
  std::vector<std::vector<unsigned int> >
    requested_ids(libMesh::n_processors());

  // We know how many of our objects live on each processor, so
  // reserve() space for requests from each.
  for (unsigned int p=0; p != libMesh::n_processors(); ++p)
    if (p != libMesh::processor_id())
      requested_ids[p].reserve(ghost_objects_from_proc[p]);

  for (it = objects_begin; it != objects_end; ++it)
    {
      DofObject *obj = *it;
      if (obj->processor_id() != DofObject::invalid_processor_id)
        requested_ids[obj->processor_id()].push_back(obj->id());
    }
#ifdef DEBUG
  for (unsigned int p=0; p != libMesh::n_processors(); ++p)
    libmesh_assert(requested_ids[p].size() == ghost_objects_from_proc[p]);
#endif

  // Next set ghost object n_comps from other processors
  for (unsigned int p=1; p != libMesh::n_processors(); ++p)
    {
      // Trade my requests with processor procup and procdown
      unsigned int procup = (libMesh::processor_id() + p) %
                             libMesh::n_processors();
      unsigned int procdown = (libMesh::n_processors() +
                               libMesh::processor_id() - p) %
                               libMesh::n_processors();
      std::vector<unsigned int> request_to_fill;
      Parallel::send_receive(procup, requested_ids[procup],
                             procdown, request_to_fill);

      // Fill those requests
      const unsigned int n_variables = this->n_variables();

      std::vector<unsigned int> ghost_data
        (request_to_fill.size() * 2 * n_variables);

      for (unsigned int i=0; i != request_to_fill.size(); ++i)
        {
          DofObject *requested = (this->*objects)(mesh, request_to_fill[i]);
          libmesh_assert(requested);
          libmesh_assert(requested->processor_id() == libMesh::processor_id());
          libmesh_assert(requested->n_vars(this->sys_number()) == n_variables);
          for (unsigned int v=0; v != n_variables; ++v)
            {
              unsigned int n_comp =
                requested->n_comp(this->sys_number(), v);
              ghost_data[i*2*n_variables+v] = n_comp;
              unsigned int first_dof = n_comp ?
                requested->dof_number(this->sys_number(), v, 0) : 0;
              libmesh_assert(first_dof != DofObject::invalid_id);
              ghost_data[i*2*n_variables+n_variables+v] = first_dof;
            }
        }

      // Trade back the results
      std::vector<unsigned int> filled_request;
      Parallel::send_receive(procdown, ghost_data,
                             procup, filled_request);

      // And copy the id changes we've now been informed of
      libmesh_assert(filled_request.size() ==
             requested_ids[procup].size() * 2 * n_variables);
      for (unsigned int i=0; i != requested_ids[procup].size(); ++i)
        {
          DofObject *requested = (this->*objects)(mesh, requested_ids[procup][i]);
          libmesh_assert(requested);
          libmesh_assert (requested->processor_id() == procup);
          for (unsigned int v=0; v != n_variables; ++v)
            {
              unsigned int n_comp = filled_request[i*2*n_variables+v];
              requested->set_n_comp(this->sys_number(), v, n_comp);
              if (n_comp)
                {
                  unsigned int first_dof = 
                    filled_request[i*2*n_variables+n_variables+v];
                  libmesh_assert(first_dof != DofObject::invalid_id);
                  requested->set_dof_number
                    (this->sys_number(), v, 0, first_dof);

                  // don't forget to add these remote dofs to the _send_list
                  for (unsigned int comp=0; comp!=n_comp; ++comp)
                    _send_list.push_back(first_dof+comp);
                }
            }
        }
    }

#ifdef DEBUG
  // Double check for invalid dofs
  for (it = objects_begin; it != objects_end; ++it)
    {
      DofObject *obj = *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
}

unsigned int DofMap::sys_number

(

)

const [inline, private]

Returns:

the number of the system we are responsible for.

Definition at line 921 of file dof_map.h.

References _sys_number.

Referenced by constrain_p_dofs(), distribute_local_dofs_node_major(), distribute_local_dofs_var_major(), dof_indices(), invalidate_dofs(), old_dof_indices(), reinit(), and set_nonlocal_dof_objects().

{
  return _sys_number;
}

bool DofMap::use_coupled_neighbor_dofs

(

const MeshBase &

mesh

)

const

Tells other library functions whether or not this problem includes coupling between dofs in neighboring cells, as can currently be specified on the command line or inferred from the use of all discontinuous variables.

Definition at line 1111 of file dof_map.C.

References FEBase::build(), libMeshEnums::DISCONTINUOUS, MeshBase::mesh_dimension(), n_variables(), libMesh::on_command_line(), and variable_type().

Referenced by compute_sparsity().

{
  // If we were asked on the command line, then we need to
  // include sensitivities between neighbor degrees of freedom
  bool implicit_neighbor_dofs = 
    libMesh::on_command_line ("--implicit_neighbor_dofs");

  // look at all the variables in this system.  If every one is
  // discontinuous then the user must be doing DG/FVM, so be nice
  // and  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;
  }
  
  return implicit_neighbor_dofs;
}

const System::Variable & DofMap::variable

(

const unsigned int

c

)

const

Returns:

the variable description object for variable c.

Definition at line 63 of file dof_map.C.

References _variables.

Referenced by distribute_local_dofs_node_major(), distribute_local_dofs_var_major(), dof_indices(), old_dof_indices(), System::ProjectVector::operator()(), reinit(), and SCALAR_dof_indices().

{
  libmesh_assert (c < _variables.size());

  return _variables[c];
}

unsigned int DofMap::variable_first_local_dof

(

const unsigned int

var

)

const [inline]

Returns the first local degree of freedom index for variable var.

Definition at line 384 of file dof_map.h.

References _var_first_local_df, and DofObject::invalid_id.

  { 
    libmesh_assert ((var+1) < _var_first_local_df.size());
    libmesh_assert (_var_first_local_df[var] != DofObject::invalid_id);
    return _var_first_local_df[var];
  }

unsigned int DofMap::variable_last_local_dof

(

const unsigned int

var

)

const [inline]

Returns (one past) the last local degree of freedom index for variable var. Analogous to the end() member function of STL containers.

Definition at line 395 of file dof_map.h.

References _var_first_local_df, and DofObject::invalid_id.

  { 
    libmesh_assert ((var+1) < _var_first_local_df.size());
    libmesh_assert (_var_first_local_df[var+1] != DofObject::invalid_id);
    return _var_first_local_df[var+1];
  }

Order DofMap::variable_order

(

const unsigned int

c

)

const

Returns:

the approximation order for variable c.

Definition at line 72 of file dof_map.C.

References _variables.

{
  libmesh_assert (c < _variables.size());

  return _variables[c].type().order;
}

const FEType & DofMap::variable_type

(

const unsigned int

c

)

const

Returns:

the finite element type for variable c.

Definition at line 81 of file dof_map.C.

References _variables.

Referenced by ExactSolution::_compute_error(), HPCoarsenTest::add_projection(), System::calculate_norm(), FEBase::coarsened_dof_values(), FEInterface::compute_constraints(), FEBase::compute_periodic_constraints(), FEBase::compute_proj_constraints(), constrain_p_dofs(), dof_indices(), JumpErrorEstimator::estimate_error(), ExactErrorEstimator::estimate_error(), System::get_info(), MeshFunction::gradient(), MeshFunction::hessian(), old_dof_indices(), System::ProjectVector::operator()(), PatchRecoveryErrorEstimator::EstimateError::operator()(), MeshFunction::operator()(), reinit(), HPCoarsenTest::select_refinement(), and use_coupled_neighbor_dofs().

{
  libmesh_assert (c < _variables.size());

  return _variables[c].type();
}

---

Friends And Related Function Documentation

friend class SparsityPattern::Build [friend]

Definition at line 870 of file dof_map.h.

---

Member Data Documentation

ReferenceCounter::Counts ReferenceCounter::_counts [static, protected, inherited]

Actually holds the data.

Definition at line 110 of file reference_counter.h.

Referenced by ReferenceCounter::get_info(), ReferenceCounter::increment_constructor_count(), and ReferenceCounter::increment_destructor_count().

DofConstraints DofMap::_dof_constraints [private]

Data structure containing DOF constraints. The ith entry is the constraint matrix row for DOF i.

Definition at line 857 of file dof_map.h.

Referenced by add_constraint_row(), add_constraints_to_send_list(), allgather_recursive_constraints(), build_constraint_matrix(), clear(), constrain_element_matrix(), constrain_element_matrix_and_vector(), constrain_p_dofs(), create_dof_constraints(), enforce_constraints_exactly(), find_connected_dofs(), is_constrained_dof(), n_constrained_dofs(), print_dof_constraints(), and process_constraints().

CouplingMatrix* DofMap::_dof_coupling

Degree of freedom coupling. If left empty each DOF couples to all others. Can be used to reduce memory requirements for sparse matrices. DOF 0 might only couple to itself, in which case dof_coupling(0,0) should be 1 and dof_coupling(0,j) = 0 for j not equal to 0.

This variable is named as though it were class private, but it is in the public interface. Also there are no public methods for accessing it... This typically means you should only use it if you know what you are doing.

Definition at line 677 of file dof_map.h.

Referenced by compute_sparsity().

std::vector<unsigned int> DofMap::_end_df [private]

Last DOF index (plus 1) on processor p.

Definition at line 806 of file dof_map.h.

Referenced by allgather_recursive_constraints(), clear(), distribute_dofs(), end_dof(), last_dof(), and n_dofs_on_processor().

std::vector<unsigned int> DofMap::_first_df [private]

First DOF index on processor p.

Definition at line 801 of file dof_map.h.

Referenced by clear(), distribute_dofs(), first_dof(), and n_dofs_on_processor().

std::vector<SparseMatrix<Number>* > DofMap::_matrices [private]

Additional matrices handled by this object. These pointers do not handle the memory, instead, System, who told DofMap about them, owns them.

Definition at line 796 of file dof_map.h.

Referenced by attach_matrix(), clear(), compute_sparsity(), and DofMap().

Threads::spin_mutex ReferenceCounter::_mutex [static, protected, inherited]

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 123 of file reference_counter.h.

unsigned int DofMap::_n_dfs [private]

Total number of degrees of freedom.

Definition at line 834 of file dof_map.h.

Referenced by clear(), distribute_dofs(), and n_dofs().

std::vector<unsigned int> DofMap::_n_nz [private]

The number of on-processor nonzeros in my portion of the global matrix.

Definition at line 823 of file dof_map.h.

Referenced by clear(), compute_sparsity(), and get_n_nz().

Threads::atomic< unsigned int > ReferenceCounter::_n_objects [static, protected, inherited]

The number of objects. Print the reference count information when the number returns to 0.

Definition at line 118 of file reference_counter.h.

Referenced by ReferenceCounter::n_objects(), ReferenceCounter::ReferenceCounter(), and ReferenceCounter::~ReferenceCounter().

unsigned int DofMap::_n_old_dfs [private]

Total number of degrees of freedom on old dof objects

Definition at line 847 of file dof_map.h.

Referenced by clear(), distribute_dofs(), and n_old_dofs().

std::vector<unsigned int> DofMap::_n_oz [private]

The number of off-processor nonzeros in my portion of the global matrix.

Definition at line 829 of file dof_map.h.

Referenced by clear(), compute_sparsity(), and get_n_oz().

unsigned int DofMap::_n_SCALAR_dofs [private]

The total number of SCALAR dofs associated to all SCALAR variables.

Definition at line 840 of file dof_map.h.

Referenced by distribute_local_dofs_node_major(), distribute_local_dofs_var_major(), n_SCALAR_dofs(), and reinit().

PeriodicBoundaries DofMap::_periodic_boundaries [private]

Data structure containing periodic boundaries. The ith entry is the constraint matrix row for boundaryid i.

Definition at line 867 of file dof_map.h.

Referenced by add_periodic_boundary(), create_dof_constraints(), and is_periodic_boundary().

std::vector<unsigned int> DofMap::_send_list [private]

A list containing all the global DOF indicies that affect the solution on my subdomain.

Definition at line 817 of file dof_map.h.

Referenced by add_constraints_to_send_list(), add_neighbors_to_send_list(), clear(), distribute_dofs(), get_send_list(), prepare_send_list(), and set_nonlocal_dof_objects().

const unsigned int DofMap::_sys_number [private]

The number of the system we manage DOFs for.

Definition at line 789 of file dof_map.h.

Referenced by sys_number().

std::vector<unsigned int> DofMap::_var_first_local_df [private]

The first local DOF index for each variable in the System.

Definition at line 811 of file dof_map.h.

Referenced by clear(), distribute_local_dofs_node_major(), distribute_local_dofs_var_major(), variable_first_local_dof(), and variable_last_local_dof().

std::vector<System::Variable> DofMap::_variables [private]

The finite element type for each variable.

Definition at line 784 of file dof_map.h.

Referenced by add_variable(), clear(), n_variables(), variable(), variable_order(), and variable_type().

---

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

• dof_map.h
• dof_map.C
• dof_map_constraints.C