libMesh::DofMap Class Reference

#include <dof_map.h>

Inheritance diagram for libMesh::DofMap:

Classes

class  AugmentSendList
 
class  AugmentSparsityPattern
 

Public Member Functions

 DofMap (const unsigned int sys_number, const ParallelObject &parent_decomp)
 
 ~DofMap ()
 
void attach_matrix (SparseMatrix< Number > &matrix)
 
bool is_attached (SparseMatrix< Number > &matrix)
 
void distribute_dofs (MeshBase &)
 
void compute_sparsity (const MeshBase &)
 
void clear_sparsity ()
 
void attach_extra_sparsity_object (DofMap::AugmentSparsityPattern &asp)
 
void attach_extra_sparsity_function (void(*func)(SparsityPattern::Graph &sparsity, std::vector< dof_id_type > &n_nz, std::vector< dof_id_type > &n_oz, void *), void *context=NULL)
 
void attach_extra_send_list_object (DofMap::AugmentSendList &asl)
 
void attach_extra_send_list_function (void(*func)(std::vector< dof_id_type > &, void *), void *context=NULL)
 
void prepare_send_list ()
 
const std::vector< dof_id_type > & get_send_list () const
 
const std::vector< dof_id_type > & get_n_nz () const
 
const std::vector< dof_id_type > & get_n_oz () const
 
void add_variable_group (const VariableGroup &var_group)
 
const VariableGroupvariable_group (const unsigned int c) const
 
const Variablevariable (const unsigned int c) const
 
Order variable_order (const unsigned int c) const
 
Order variable_group_order (const unsigned int vg) const
 
const FETypevariable_type (const unsigned int c) const
 
const FETypevariable_group_type (const unsigned int vg) const
 
unsigned int n_variable_groups () const
 
unsigned int n_variables () const
 
bool has_blocked_representation () const
 
unsigned int block_size () const
 
dof_id_type n_dofs () const
 
dof_id_type n_SCALAR_dofs () const
 
dof_id_type n_local_dofs () const
 
dof_id_type n_dofs_on_processor (const processor_id_type proc) const
 
dof_id_type first_dof (const processor_id_type proc) const
 
dof_id_type first_dof () const
 
dof_id_type first_old_dof (const processor_id_type proc) const
 
dof_id_type first_old_dof () const
 
dof_id_type last_dof (const processor_id_type proc) const
 
dof_id_type last_dof () const
 
dof_id_type end_dof (const processor_id_type proc) const
 
dof_id_type end_dof () const
 
dof_id_type end_old_dof (const processor_id_type proc) const
 
dof_id_type end_old_dof () const
 
void dof_indices (const Elem *const elem, std::vector< dof_id_type > &di, const unsigned int vn=libMesh::invalid_uint) const
 
void SCALAR_dof_indices (std::vector< dof_id_type > &di, const unsigned int vn, const bool old_dofs=false) const
 
bool all_semilocal_indices (const std::vector< dof_id_type > &dof_indices) const
 
bool use_coupled_neighbor_dofs (const MeshBase &mesh) const
 
void extract_local_vector (const NumericVector< Number > &Ug, const std::vector< dof_id_type > &dof_indices, DenseVectorBase< Number > &Ue) const
 
dof_id_type n_constrained_dofs () const
 
dof_id_type n_local_constrained_dofs () const
 
dof_id_type n_constrained_nodes () const
 
void create_dof_constraints (const MeshBase &, Real time=0)
 
void allgather_recursive_constraints (MeshBase &)
 
void scatter_constraints (MeshBase &)
 
void process_constraints (MeshBase &)
 
void add_constraint_row (const dof_id_type dof_number, const DofConstraintRow &constraint_row, const Number constraint_rhs, const bool forbid_constraint_overwrite)
 
void add_adjoint_constraint_row (const unsigned int qoi_index, const dof_id_type dof_number, const DofConstraintRow &constraint_row, const Number constraint_rhs, const bool forbid_constraint_overwrite)
 
void add_constraint_row (const dof_id_type dof_number, const DofConstraintRow &constraint_row, const bool forbid_constraint_overwrite=true)
 
DofConstraints::const_iterator constraint_rows_begin () const
 
DofConstraints::const_iterator constraint_rows_end () const
 
NodeConstraints::const_iterator node_constraint_rows_begin () const
 
NodeConstraints::const_iterator node_constraint_rows_end () const
 
bool is_constrained_dof (const dof_id_type dof) const
 
bool has_heterogenous_adjoint_constraint (const unsigned int qoi_num, const dof_id_type dof) const
 
bool is_constrained_node (const Node *node) const
 
void print_dof_constraints (std::ostream &os=libMesh::out, bool print_nonlocal=false) const
 
std::string get_local_constraints (bool print_nonlocal=false) const
 
std::pair< Real, Realmax_constraint_error (const System &system, NumericVector< Number > *v=NULL) const
 
void constrain_element_matrix (DenseMatrix< Number > &matrix, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
 
void constrain_element_matrix (DenseMatrix< Number > &matrix, std::vector< dof_id_type > &row_dofs, std::vector< dof_id_type > &col_dofs, bool asymmetric_constraint_rows=true) const
 
void constrain_element_vector (DenseVector< Number > &rhs, std::vector< dof_id_type > &dofs, bool asymmetric_constraint_rows=true) const
 
void constrain_element_matrix_and_vector (DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
 
void heterogenously_constrain_element_matrix_and_vector (DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true, int qoi_index=-1) const
 
void heterogenously_constrain_element_vector (const DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true, int qoi_index=-1) const
 
void constrain_element_dyad_matrix (DenseVector< Number > &v, DenseVector< Number > &w, std::vector< dof_id_type > &row_dofs, bool asymmetric_constraint_rows=true) const
 
void constrain_nothing (std::vector< dof_id_type > &dofs) const
 
void enforce_constraints_exactly (const System &system, NumericVector< Number > *v=NULL, bool homogeneous=false) const
 
void enforce_adjoint_constraints_exactly (NumericVector< Number > &v, unsigned int q) const
 
void add_periodic_boundary (const PeriodicBoundaryBase &periodic_boundary)
 
void add_periodic_boundary (const PeriodicBoundaryBase &boundary, const PeriodicBoundaryBase &inverse_boundary)
 
bool is_periodic_boundary (const boundary_id_type boundaryid) const
 
PeriodicBoundariesget_periodic_boundaries ()
 
void add_dirichlet_boundary (const DirichletBoundary &dirichlet_boundary)
 
void add_adjoint_dirichlet_boundary (const DirichletBoundary &dirichlet_boundary, unsigned int q)
 
void remove_dirichlet_boundary (const DirichletBoundary &dirichlet_boundary)
 
void remove_adjoint_dirichlet_boundary (const DirichletBoundary &dirichlet_boundary, unsigned int q)
 
const DirichletBoundariesget_dirichlet_boundaries () const
 
DirichletBoundariesget_dirichlet_boundaries ()
 
bool has_adjoint_dirichlet_boundaries (unsigned int q) const
 
const DirichletBoundariesget_adjoint_dirichlet_boundaries (unsigned int q) const
 
DirichletBoundariesget_adjoint_dirichlet_boundaries (unsigned int q)
 
void old_dof_indices (const Elem *const elem, std::vector< dof_id_type > &di, const unsigned int vn=libMesh::invalid_uint) const
 
dof_id_type 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 ()
 
void print_info (std::ostream &os=libMesh::out) const
 
std::string get_info () const
 
unsigned int sys_number () const
 
const Parallel::Communicatorcomm () const
 
processor_id_type n_processors () const
 
processor_id_type processor_id () const
 

Static Public Member Functions

static std::string get_info ()
 
static void print_info (std::ostream &out=libMesh::out)
 
static unsigned int n_objects ()
 
static void enable_print_counter_info ()
 
static void disable_print_counter_info ()
 

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)
 

Protected Attributes

const Parallel::Communicator_communicator
 

Static Protected Attributes

static Counts _counts
 
static Threads::atomic
< unsigned int > 
_n_objects
 
static Threads::spin_mutex _mutex
 
static bool _enable_print_counter = true
 

Private Types

typedef DofObject *(DofMap::* dofobject_accessor )(MeshBase &mesh, dof_id_type i) const
 

Private Member Functions

AutoPtr< SparsityPattern::Buildbuild_sparsity (const MeshBase &mesh) const
 
void invalidate_dofs (MeshBase &mesh) const
 
DofObjectnode_ptr (MeshBase &mesh, dof_id_type i) const
 
DofObjectelem_ptr (MeshBase &mesh, dof_id_type 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 (dof_id_type &next_free_dof, MeshBase &mesh)
 
void distribute_local_dofs_node_major (dof_id_type &next_free_dof, MeshBase &mesh)
 
void add_neighbors_to_send_list (MeshBase &mesh)
 
void build_constraint_matrix (DenseMatrix< Number > &C, std::vector< dof_id_type > &elem_dofs, const bool called_recursively=false) const
 
void build_constraint_matrix_and_vector (DenseMatrix< Number > &C, DenseVector< Number > &H, std::vector< dof_id_type > &elem_dofs, int qoi_index=-1, const bool called_recursively=false) const
 
void find_connected_dofs (std::vector< dof_id_type > &elem_dofs) const
 
void find_connected_dof_objects (std::vector< const DofObject * > &objs) const
 
void add_constraints_to_send_list ()
 

Private Attributes

std::vector< Variable_variables
 
std::vector< VariableGroup_variable_groups
 
const unsigned int _sys_number
 
std::vector< SparseMatrix
< Number > * > 
_matrices
 
std::vector< dof_id_type_first_df
 
std::vector< dof_id_type_end_df
 
std::vector< dof_id_type_send_list
 
AugmentSparsityPattern_augment_sparsity_pattern
 
void(* _extra_sparsity_function )(SparsityPattern::Graph &, std::vector< dof_id_type > &n_nz, std::vector< dof_id_type > &n_oz, void *)
 
void * _extra_sparsity_context
 
AugmentSendList_augment_send_list
 
void(* _extra_send_list_function )(std::vector< dof_id_type > &, void *)
 
void * _extra_send_list_context
 
bool need_full_sparsity_pattern
 
AutoPtr< SparsityPattern::Build_sp
 
std::vector< dof_id_type > * _n_nz
 
std::vector< dof_id_type > * _n_oz
 
dof_id_type _n_dfs
 
dof_id_type _n_SCALAR_dofs
 
dof_id_type _n_old_dfs
 
std::vector< dof_id_type_first_old_df
 
std::vector< dof_id_type_end_old_df
 
DofConstraints _dof_constraints
 
DofConstraintValueMap _primal_constraint_values
 
AdjointDofConstraintValues _adjoint_constraint_values
 
NodeConstraints _node_constraints
 
PeriodicBoundaries_periodic_boundaries
 
DirichletBoundaries_dirichlet_boundaries
 
std::vector
< DirichletBoundaries * > 
_adjoint_dirichlet_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 167 of file dof_map.h.

Member Typedef Documentation

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > libMesh::ReferenceCounter::Counts
protectedinherited

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

Definition at line 113 of file reference_counter.h.

typedef DofObject*(DofMap::* libMesh::DofMap::dofobject_accessor)(MeshBase &mesh, dof_id_type i) const
private

A member function type like node_ptr or elem_ptr

Definition at line 1071 of file dof_map.h.

Constructor & Destructor Documentation

libMesh::DofMap::DofMap ( const unsigned int  sys_number,
const ParallelObject parent_decomp 
)
explicit

Constructor. Requires the number of the system for which we will be numbering degrees of freedom & the parent object we are contained in, which defines our communication space.

Definition at line 177 of file dof_map.C.

References _adjoint_dirichlet_boundaries, _dirichlet_boundaries, _periodic_boundaries, and clear().

178 {
179  this->clear();
180 #ifdef LIBMESH_ENABLE_PERIODIC
181  delete _periodic_boundaries;
182 #endif
183 #ifdef LIBMESH_ENABLE_DIRICHLET
184  delete _dirichlet_boundaries;
185  for (unsigned int q = 0; q != _adjoint_dirichlet_boundaries.size(); ++q)
187 #endif
188 }
libMesh::DofMap::~DofMap ( )

Destructor.

Member Function Documentation

void libMesh::DofMap::add_adjoint_constraint_row ( const unsigned int  qoi_index,
const dof_id_type  dof_number,
const DofConstraintRow constraint_row,
const Number  constraint_rhs,
const bool  forbid_constraint_overwrite 
)

Adds a copy of the user-defined row to the constraint matrix, using an inhomogeneous right-hand-side for the adjoint constraint equation.

forbid_constraint_overwrite here only tests for overwriting the rhs. This method should only be used when an equivalent constraint (with a potentially different rhs) already exists for the primal problem.

Definition at line 1116 of file dof_map_constraints.C.

References libMesh::err.

1121 {
1122  // Optionally allow the user to overwrite constraints. Defaults to false.
1123  if (forbid_constraint_overwrite)
1124  {
1125  if (!this->is_constrained_dof(dof_number))
1126  {
1127  libMesh::err << "ERROR: DOF " << dof_number << " has no corresponding primal constraint!"
1128  << std::endl;
1129  libmesh_error();
1130  }
1131 #ifndef NDEBUG
1132  // No way to do this without a non-normalized tolerance?
1133  /*
1134  // If the user passed in more than just the rhs, let's check the
1135  // coefficients for consistency
1136  if (!constraint_row.empty())
1137  {
1138  DofConstraintRow row = _dof_constraints[dof_number];
1139  for (DofConstraintRow::const_iterator pos=row.begin();
1140  pos != row.end(); ++pos)
1141  {
1142  libmesh_assert(constraint_row.find(pos->first)->second
1143  == pos->second);
1144  }
1145  }
1146 
1147  if (_adjoint_constraint_values[qoi_index].find(dof_number) !=
1148  _adjoint_constraint_values[qoi_index].end())
1149  libmesh_assert_equal_to
1150  (_adjoint_constraint_values[qoi_index][dof_number],
1151  constraint_rhs);
1152  */
1153 #endif
1154  }
1155 
1156  // Creates the map of rhs values if it doesn't already exist; then
1157  // adds the current value to that map
1158  _adjoint_constraint_values[qoi_index].insert(std::make_pair(dof_number, constraint_rhs));
1159 }
void libMesh::DofMap::add_adjoint_dirichlet_boundary ( const DirichletBoundary dirichlet_boundary,
unsigned int  q 
)

Adds a copy of the specified Dirichlet boundary to the system, corresponding to the adjoint problem defined by Quantity of Interest q.

Definition at line 3562 of file dof_map_constraints.C.

3564 {
3565  std::size_t old_size = _adjoint_dirichlet_boundaries.size();
3566  for (std::size_t i = old_size; i <= qoi_index; ++i)
3568 
3569  _adjoint_dirichlet_boundaries[qoi_index]->push_back(new DirichletBoundary(dirichlet_boundary));
3570 }
void libMesh::DofMap::add_constraint_row ( const dof_id_type  dof_number,
const DofConstraintRow constraint_row,
const Number  constraint_rhs,
const bool  forbid_constraint_overwrite 
)

Adds a copy of the user-defined row to the constraint matrix, using an inhomogeneous right-hand-side for the constraint equation.

Definition at line 1097 of file dof_map_constraints.C.

References libMesh::err.

Referenced by add_constraint_row().

1101 {
1102  // Optionally allow the user to overwrite constraints. Defaults to false.
1103  if (forbid_constraint_overwrite)
1104  if (this->is_constrained_dof(dof_number))
1105  {
1106  libMesh::err << "ERROR: DOF " << dof_number << " was already constrained!"
1107  << std::endl;
1108  libmesh_error();
1109  }
1110 
1111  _dof_constraints.insert(std::make_pair(dof_number, constraint_row));
1112  _primal_constraint_values.insert(std::make_pair(dof_number, constraint_rhs));
1113 }
void libMesh::DofMap::add_constraint_row ( const dof_id_type  dof_number,
const DofConstraintRow constraint_row,
const bool  forbid_constraint_overwrite = true 
)
inline

Adds a copy of the user-defined row to the constraint matrix, using a homogeneous right-hand-side for the constraint equation. By default, produces an error if the DOF was already constrained.

Definition at line 638 of file dof_map.h.

References add_constraint_row().

641  { add_constraint_row(dof_number, constraint_row, 0., forbid_constraint_overwrite); }
void libMesh::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 3440 of file dof_map_constraints.C.

References libMesh::comm, and libMesh::n_processors().

3441 {
3442  // This function must be run on all processors at once
3443  parallel_object_only();
3444 
3445  // Return immediately if there's nothing to gather
3446  if (this->n_processors() == 1)
3447  return;
3448 
3449  // We might get to return immediately if none of the processors
3450  // found any constraints
3451  unsigned int has_constraints = !_dof_constraints.empty();
3452  this->comm().max(has_constraints);
3453  if (!has_constraints)
3454  return;
3455 
3456  for (DofConstraints::iterator i = _dof_constraints.begin();
3457  i != _dof_constraints.end(); ++i)
3458  {
3459  dof_id_type constrained_dof = i->first;
3460 
3461  // We only need the dependencies of our own constrained dofs
3462  if (constrained_dof < this->first_dof() ||
3463  constrained_dof >= this->end_dof())
3464  continue;
3465 
3466  DofConstraintRow& constraint_row = i->second;
3467  for (DofConstraintRow::const_iterator
3468  j=constraint_row.begin(); j != constraint_row.end();
3469  ++j)
3470  {
3471  dof_id_type constraint_dependency = j->first;
3472 
3473  // No point in adding one of our own dofs to the send_list
3474  if (constraint_dependency >= this->first_dof() &&
3475  constraint_dependency < this->end_dof())
3476  continue;
3477 
3478  _send_list.push_back(constraint_dependency);
3479  }
3480  }
3481 }
void libMesh::DofMap::add_dirichlet_boundary ( const DirichletBoundary dirichlet_boundary)

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

Definition at line 3555 of file dof_map_constraints.C.

3556 {
3557  _dirichlet_boundaries->push_back(new DirichletBoundary(dirichlet_boundary));
3558 }
void libMesh::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 1293 of file dof_map.C.

References _send_list, libMesh::Elem::active(), libMesh::MeshBase::active_elements_begin(), libMesh::MeshBase::active_elements_end(), libMesh::Elem::active_family_tree_by_neighbor(), libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), dof_indices(), libMesh::DofObject::dof_number(), end_dof(), first_dof(), libMesh::Elem::get_node(), libMesh::MeshBase::max_node_id(), libMesh::DofObject::n_comp(), libMesh::Elem::n_neighbors(), libMesh::Elem::n_nodes(), libMesh::n_vars, libMesh::DofObject::n_vars(), libMesh::Elem::neighbor(), libMesh::Elem::node(), libMesh::ParallelObject::processor_id(), libMesh::DofObject::processor_id(), libMesh::START_LOG(), libMesh::STOP_LOG(), and sys_number().

Referenced by distribute_dofs().

1294 {
1295  START_LOG("add_neighbors_to_send_list()", "DofMap");
1296 
1297  const unsigned int sys_num = this->sys_number();
1298 
1299  //-------------------------------------------------------------------------
1300  // We need to add the DOFs from elements that live on neighboring processors
1301  // that are neighbors of the elements on the local processor
1302  //-------------------------------------------------------------------------
1303 
1304  MeshBase::const_element_iterator local_elem_it
1305  = mesh.active_local_elements_begin();
1306  const MeshBase::const_element_iterator local_elem_end
1307  = mesh.active_local_elements_end();
1308 
1309  std::vector<bool> node_on_processor(mesh.max_node_id(), false);
1310  std::vector<dof_id_type> di;
1311  std::vector<const Elem *> family;
1312 
1313  // Loop over the active local elements, adding all active elements
1314  // that neighbor an active local element to the send list.
1315  for ( ; local_elem_it != local_elem_end; ++local_elem_it)
1316  {
1317  const Elem* elem = *local_elem_it;
1318 
1319  for (unsigned int n=0; n!=elem->n_nodes(); n++)
1320  {
1321  // Flag all the nodes of active local elements as seen, so
1322  // we can add nodal neighbor dofs to the send_list later.
1323  node_on_processor[elem->node(n)] = true;
1324 
1325  // Add all remote dofs on these nodes to the send_list.
1326  // This is necessary in case those dofs are *not* also dofs
1327  // on neighbors; e.g. in the case of a HIERARCHIC's local
1328  // side which is only a vertex on the neighbor that owns it.
1329  const Node* node = elem->get_node(n);
1330  const unsigned n_vars = node->n_vars(sys_num);
1331  for (unsigned int v=0; v != n_vars; ++v)
1332  {
1333  const unsigned int n_comp = node->n_comp(sys_num, v);
1334  for (unsigned int c=0; c != n_comp; ++c)
1335  {
1336  const dof_id_type dof_index = node->dof_number(sys_num, v, c);
1337  if (dof_index < this->first_dof() || dof_index >= this->end_dof())
1338  {
1339  _send_list.push_back(dof_index);
1340  // libmesh_here();
1341  // libMesh::out << "sys_num,v,c,dof_index="
1342  // << sys_num << ", "
1343  // << v << ", "
1344  // << c << ", "
1345  // << dof_index << '\n';
1346  // node->debug_buffer();
1347  }
1348  }
1349  }
1350  }
1351 
1352  // Loop over the neighbors of those elements
1353  for (unsigned int s=0; s<elem->n_neighbors(); s++)
1354  if (elem->neighbor(s) != NULL)
1355  {
1356  family.clear();
1357 
1358  // Find all the active elements that neighbor elem
1359 #ifdef LIBMESH_ENABLE_AMR
1360  if (!elem->neighbor(s)->active())
1361  elem->neighbor(s)->active_family_tree_by_neighbor(family, elem);
1362  else
1363 #endif
1364  family.push_back(elem->neighbor(s));
1365 
1366  for (dof_id_type i=0; i!=family.size(); ++i)
1367  // If the neighbor lives on a different processor
1368  if (family[i]->processor_id() != this->processor_id())
1369  {
1370  // Get the DOF indices for this neighboring element
1371  this->dof_indices (family[i], di);
1372 
1373  // Insert the remote DOF indices into the send list
1374  for (std::size_t j=0; j != di.size(); ++j)
1375  if (di[j] < this->first_dof() ||
1376  di[j] >= this->end_dof())
1377  _send_list.push_back(di[j]);
1378  }
1379  }
1380  }
1381 
1382  // Now loop over all non_local active elements and add any missing
1383  // nodal-only neighbors. This will also get any dofs from nonlocal
1384  // nodes on local elements, because every nonlocal node exists on a
1385  // nonlocal nodal neighbor element.
1386  MeshBase::const_element_iterator elem_it
1387  = mesh.active_elements_begin();
1388  const MeshBase::const_element_iterator elem_end
1389  = mesh.active_elements_end();
1390 
1391  for ( ; elem_it != elem_end; ++elem_it)
1392  {
1393  const Elem* elem = *elem_it;
1394 
1395  // If this is one of our elements, we've already added it
1396  if (elem->processor_id() == this->processor_id())
1397  continue;
1398 
1399  // Do we need to add the element DOFs?
1400  bool add_elem_dofs = false;
1401 
1402  // Check all the nodes of the element to see if it
1403  // shares a node with us
1404  for (unsigned int n=0; n!=elem->n_nodes(); n++)
1405  if (node_on_processor[elem->node(n)])
1406  add_elem_dofs = true;
1407 
1408  // Add the element degrees of freedom if it shares at
1409  // least one node.
1410  if (add_elem_dofs)
1411  {
1412  // Get the DOF indices for this neighboring element
1413  this->dof_indices (elem, di);
1414 
1415  // Insert the remote DOF indices into the send list
1416  for (std::size_t j=0; j != di.size(); ++j)
1417  if (di[j] < this->first_dof() ||
1418  di[j] >= this->end_dof())
1419  _send_list.push_back(di[j]);
1420  }
1421  }
1422 
1423  STOP_LOG("add_neighbors_to_send_list()", "DofMap");
1424 }
void libMesh::DofMap::add_periodic_boundary ( const PeriodicBoundaryBase periodic_boundary)

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

Definition at line 3660 of file dof_map_constraints.C.

References libMesh::PeriodicBoundaryBase::clone(), libMesh::libmesh_assert(), libMesh::PeriodicBoundaryBase::merge(), libMesh::PeriodicBoundaryBase::myboundary, and libMesh::PeriodicBoundaryBase::pairedboundary.

3661 {
3662  // See if we already have a periodic boundary associated myboundary...
3663  PeriodicBoundaryBase* existing_boundary = _periodic_boundaries->boundary(periodic_boundary.myboundary);
3664 
3665  if ( existing_boundary == NULL )
3666  {
3667  // ...if not, clone the input (and its inverse) and add them to the PeriodicBoundaries object
3668  PeriodicBoundaryBase* boundary = periodic_boundary.clone().release();
3669  PeriodicBoundaryBase* inverse_boundary = periodic_boundary.clone(PeriodicBoundaryBase::INVERSE).release();
3670 
3671  // _periodic_boundaries takes ownership of the pointers
3672  _periodic_boundaries->insert(std::make_pair(boundary->myboundary, boundary));
3673  _periodic_boundaries->insert(std::make_pair(inverse_boundary->myboundary, inverse_boundary));
3674  }
3675  else
3676  {
3677  // ...otherwise, merge this object's variable IDs with the existing boundary object's.
3678  existing_boundary->merge(periodic_boundary);
3679 
3680  // Do the same merging process for the inverse boundary. Note: the inverse better already exist!
3681  PeriodicBoundaryBase* inverse_boundary = _periodic_boundaries->boundary(periodic_boundary.pairedboundary);
3682  libmesh_assert(inverse_boundary);
3683  inverse_boundary->merge(periodic_boundary);
3684  }
3685 }
void libMesh::DofMap::add_periodic_boundary ( const PeriodicBoundaryBase boundary,
const PeriodicBoundaryBase inverse_boundary 
)

Add a periodic boundary pair

Parameters
boundary- primary boundary
inverse_boundary- inverse boundary

Definition at line 3690 of file dof_map_constraints.C.

References libMesh::PeriodicBoundaryBase::clone(), libMesh::PeriodicBoundaryBase::myboundary, and libMesh::PeriodicBoundaryBase::pairedboundary.

3691 {
3692  libmesh_assert_equal_to (boundary.myboundary, inverse_boundary.pairedboundary);
3693  libmesh_assert_equal_to (boundary.pairedboundary, inverse_boundary.myboundary);
3694 
3695  // Allocate copies on the heap. The _periodic_boundaries object will manage this memory.
3696  // Note: this also means that the copy constructor for the PeriodicBoundary (or user class
3697  // derived therefrom) must be implemented!
3698  // PeriodicBoundary* p_boundary = new PeriodicBoundary(boundary);
3699  // PeriodicBoundary* p_inverse_boundary = new PeriodicBoundary(inverse_boundary);
3700 
3701  // We can't use normal copy construction since this leads to slicing with derived classes.
3702  // Use clone() "virtual constructor" instead. But, this *requires* user to override the clone()
3703  // method. Note also that clone() allocates memory. In this case, the _periodic_boundaries object
3704  // takes responsibility for cleanup.
3705  PeriodicBoundaryBase* p_boundary = boundary.clone().release();
3706  PeriodicBoundaryBase* p_inverse_boundary = inverse_boundary.clone().release();
3707 
3708  // Add the periodic boundary and its inverse to the PeriodicBoundaries data structure. The
3709  // PeriodicBoundaries data structure takes ownership of the pointers.
3710  _periodic_boundaries->insert(std::make_pair(p_boundary->myboundary, p_boundary));
3711  _periodic_boundaries->insert(std::make_pair(p_inverse_boundary->myboundary, p_inverse_boundary));
3712 }
void libMesh::DofMap::add_variable_group ( const VariableGroup var_group)

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

Definition at line 213 of file dof_map.C.

References _variable_groups, _variables, and libMesh::VariableGroup::n_variables().

214 {
215  _variable_groups.push_back(var_group);
216 
217  VariableGroup &new_var_group = _variable_groups.back();
218 
219  for (unsigned int var=0; var<new_var_group.n_variables(); var++)
220  _variables.push_back (new_var_group(var));
221 }
bool libMesh::DofMap::all_semilocal_indices ( const std::vector< dof_id_type > &  dof_indices) const

Returns true iff all degree of freedom indices in dof_indices are either local indices or in the send_list. Note that this is an O(logN) operation, not O(1); we don't cache enough information for O(1) right now.

Definition at line 1863 of file dof_map.C.

References _send_list, end_dof(), and first_dof().

1864 {
1865  // We're all semilocal unless we find a counterexample
1866  for (std::size_t i=0; i != dof_indices_in.size(); ++i)
1867  {
1868  const dof_id_type di = dof_indices_in[i];
1869  // If it's not in the local indices
1870  if (di < this->first_dof() ||
1871  di >= this->end_dof())
1872  {
1873  // and if it's not in the ghost indices, then we're not
1874  // semilocal
1875  if (!std::binary_search(_send_list.begin(), _send_list.end(), di))
1876  return false;
1877  }
1878  }
1879 
1880  return true;
1881 }
void libMesh::DofMap::allgather_recursive_constraints ( MeshBase mesh)

Gathers constraint equation dependencies from other processors

Definition at line 2400 of file dof_map_constraints.C.

References libMesh::MeshBase::active_not_local_elements_begin(), libMesh::MeshBase::active_not_local_elements_end(), libMesh::comm, end, libMesh::Elem::get_node(), libMesh::DofObject::id(), libMesh::libmesh_assert(), mesh, libMesh::Elem::n_nodes(), libMesh::n_processors(), libMesh::Elem::node(), libMesh::MeshBase::node_ptr(), libMesh::processor_id(), libMesh::DofObject::processor_id(), and libMesh::TypeVector< T >::size().

2401 {
2402  // This function must be run on all processors at once
2403  parallel_object_only();
2404 
2405  // Return immediately if there's nothing to gather
2406  if (this->n_processors() == 1)
2407  return;
2408 
2409  // We might get to return immediately if none of the processors
2410  // found any constraints
2411  unsigned int has_constraints = !_dof_constraints.empty()
2412 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2413  || !_node_constraints.empty()
2414 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2415  ;
2416  this->comm().max(has_constraints);
2417  if (!has_constraints)
2418  return;
2419 
2420  // We might have calculated constraints for constrained dofs
2421  // which have support on other processors.
2422  // Push these out first.
2423  {
2424  std::vector<std::set<dof_id_type> > pushed_ids(this->n_processors());
2425 
2426 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2427  std::vector<std::set<dof_id_type> > pushed_node_ids(this->n_processors());
2428 #endif
2429 
2431  foreign_elem_it = mesh.active_not_local_elements_begin(),
2432  foreign_elem_end = mesh.active_not_local_elements_end();
2433 
2434  // Collect the constraints to push to each processor
2435  for (; foreign_elem_it != foreign_elem_end; ++foreign_elem_it)
2436  {
2437  Elem *elem = *foreign_elem_it;
2438 
2439  std::vector<dof_id_type> my_dof_indices;
2440  this->dof_indices (elem, my_dof_indices);
2441 
2442  for (unsigned int i=0; i != my_dof_indices.size(); ++i)
2443  if (this->is_constrained_dof(my_dof_indices[i]))
2444  pushed_ids[elem->processor_id()].insert(my_dof_indices[i]);
2445 
2446 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2447  for (unsigned int n=0; n != elem->n_nodes(); ++n)
2448  if (this->is_constrained_node(elem->get_node(n)))
2449  pushed_node_ids[elem->processor_id()].insert(elem->node(n));
2450 #endif
2451  }
2452 
2453  // Now trade constraint rows
2454  for (processor_id_type p = 0; p != this->n_processors(); ++p)
2455  {
2456  // Push to processor procup while receiving from procdown
2457  processor_id_type procup = (this->processor_id() + p) %
2458  this->n_processors();
2459  processor_id_type procdown = (this->n_processors() +
2460  this->processor_id() - p) %
2461  this->n_processors();
2462 
2463  // Pack the dof constraint rows and rhs's to push to procup
2464  const std::size_t pushed_ids_size = pushed_ids[procup].size();
2465  std::vector<std::vector<dof_id_type> > pushed_keys(pushed_ids_size);
2466  std::vector<std::vector<Real> > pushed_vals(pushed_ids_size);
2467  std::vector<Number> pushed_rhss(pushed_ids_size);
2468  std::set<dof_id_type>::const_iterator it = pushed_ids[procup].begin();
2469  for (std::size_t i = 0; it != pushed_ids[procup].end();
2470  ++i, ++it)
2471  {
2472  const dof_id_type pushed_id = *it;
2473  DofConstraintRow &row = _dof_constraints[pushed_id];
2474  std::size_t row_size = row.size();
2475  pushed_keys[i].reserve(row_size);
2476  pushed_vals[i].reserve(row_size);
2477  for (DofConstraintRow::iterator j = row.begin();
2478  j != row.end(); ++j)
2479  {
2480  pushed_keys[i].push_back(j->first);
2481  pushed_vals[i].push_back(j->second);
2482  }
2483 
2484  DofConstraintValueMap::const_iterator rhsit =
2485  _primal_constraint_values.find(pushed_id);
2486  pushed_rhss[i] =
2487  (rhsit == _primal_constraint_values.end()) ?
2488  0 : rhsit->second;
2489  }
2490 
2491 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2492  // Pack the node constraint rows to push to procup
2493  const std::size_t pushed_nodes_size = pushed_node_ids[procup].size();
2494  std::vector<std::vector<dof_id_type> > pushed_node_keys(pushed_nodes_size);
2495  std::vector<std::vector<Real> > pushed_node_vals(pushed_nodes_size);
2496  std::vector<Point> pushed_node_offsets(pushed_nodes_size);
2497  std::set<dof_id_type>::const_iterator node_it = pushed_node_ids[procup].begin();
2498  for (std::size_t i = 0; node_it != pushed_node_ids[procup].end();
2499  ++i, ++node_it)
2500  {
2501  const Node *node = mesh.node_ptr(*node_it);
2502  NodeConstraintRow &row = _node_constraints[node].first;
2503  std::size_t row_size = row.size();
2504  pushed_node_keys[i].reserve(row_size);
2505  pushed_node_vals[i].reserve(row_size);
2506  for (NodeConstraintRow::iterator j = row.begin();
2507  j != row.end(); ++j)
2508  {
2509  pushed_node_keys[i].push_back(j->first->id());
2510  pushed_node_vals[i].push_back(j->second);
2511  }
2512  pushed_node_offsets[i] = _node_constraints[node].second;
2513  }
2514 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2515 
2516  // Trade pushed dof constraint rows
2517  std::vector<dof_id_type> pushed_ids_from_me
2518  (pushed_ids[procup].begin(), pushed_ids[procup].end());
2519  std::vector<dof_id_type> pushed_ids_to_me;
2520  std::vector<std::vector<dof_id_type> > pushed_keys_to_me;
2521  std::vector<std::vector<Real> > pushed_vals_to_me;
2522  std::vector<Number> pushed_rhss_to_me;
2523  this->comm().send_receive(procup, pushed_ids_from_me,
2524  procdown, pushed_ids_to_me);
2525  this->comm().send_receive(procup, pushed_keys,
2526  procdown, pushed_keys_to_me);
2527  this->comm().send_receive(procup, pushed_vals,
2528  procdown, pushed_vals_to_me);
2529  this->comm().send_receive(procup, pushed_rhss,
2530  procdown, pushed_rhss_to_me);
2531  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_keys_to_me.size());
2532  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_vals_to_me.size());
2533  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_rhss_to_me.size());
2534 
2535 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2536  // Trade pushed node constraint rows
2537  std::vector<dof_id_type> pushed_node_ids_from_me
2538  (pushed_node_ids[procup].begin(), pushed_node_ids[procup].end());
2539  std::vector<dof_id_type> pushed_node_ids_to_me;
2540  std::vector<std::vector<dof_id_type> > pushed_node_keys_to_me;
2541  std::vector<std::vector<Real> > pushed_node_vals_to_me;
2542  std::vector<Point> pushed_node_offsets_to_me;
2543  this->comm().send_receive(procup, pushed_node_ids_from_me,
2544  procdown, pushed_node_ids_to_me);
2545  this->comm().send_receive(procup, pushed_node_keys,
2546  procdown, pushed_node_keys_to_me);
2547  this->comm().send_receive(procup, pushed_node_vals,
2548  procdown, pushed_node_vals_to_me);
2549  this->comm().send_receive(procup, pushed_node_offsets,
2550  procdown, pushed_node_offsets_to_me);
2551 
2552  // Note that we aren't doing a send_receive on the Nodes
2553  // themselves. At this point we should only be pushing out
2554  // "raw" constraints, and there should be no
2555  // constrained-by-constrained-by-etc. situations that could
2556  // involve non-semilocal nodes.
2557 
2558  libmesh_assert_equal_to (pushed_node_ids_to_me.size(), pushed_node_keys_to_me.size());
2559  libmesh_assert_equal_to (pushed_node_ids_to_me.size(), pushed_node_vals_to_me.size());
2560  libmesh_assert_equal_to (pushed_node_ids_to_me.size(), pushed_node_offsets_to_me.size());
2561 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2562 
2563  // Add the dof constraints that I've been sent
2564  for (std::size_t i = 0; i != pushed_ids_to_me.size(); ++i)
2565  {
2566  libmesh_assert_equal_to (pushed_keys_to_me[i].size(), pushed_vals_to_me[i].size());
2567 
2568  dof_id_type constrained = pushed_ids_to_me[i];
2569 
2570  // If we don't already have a constraint for this dof,
2571  // add the one we were sent
2572  if (!this->is_constrained_dof(constrained))
2573  {
2574  DofConstraintRow &row = _dof_constraints[constrained];
2575  for (std::size_t j = 0; j != pushed_keys_to_me[i].size(); ++j)
2576  {
2577  row[pushed_keys_to_me[i][j]] = pushed_vals_to_me[i][j];
2578  }
2579  if (pushed_rhss_to_me[i] != Number(0))
2580  _primal_constraint_values[constrained] = pushed_rhss_to_me[i];
2581  else
2582  _primal_constraint_values.erase(constrained);
2583  }
2584  }
2585 
2586 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2587  // Add the node constraints that I've been sent
2588  for (std::size_t i = 0; i != pushed_node_ids_to_me.size(); ++i)
2589  {
2590  libmesh_assert_equal_to (pushed_node_keys_to_me[i].size(), pushed_node_vals_to_me[i].size());
2591 
2592  dof_id_type constrained_id = pushed_node_ids_to_me[i];
2593 
2594  // If we don't already have a constraint for this node,
2595  // add the one we were sent
2596  const Node *constrained = mesh.node_ptr(constrained_id);
2597  if (!this->is_constrained_node(constrained))
2598  {
2599  NodeConstraintRow &row = _node_constraints[constrained].first;
2600  for (std::size_t j = 0; j != pushed_node_keys_to_me[i].size(); ++j)
2601  {
2602  const Node *key_node = mesh.node_ptr(pushed_node_keys_to_me[i][j]);
2603  libmesh_assert(key_node);
2604  row[key_node] = pushed_node_vals_to_me[i][j];
2605  }
2606  _node_constraints[constrained].second = pushed_node_offsets_to_me[i];
2607  }
2608  }
2609 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2610  }
2611  }
2612 
2613  // Now start checking for any other constraints we need
2614  // to know about, requesting them recursively.
2615 
2616  // Create sets containing the DOFs and nodes we already depend on
2617  typedef std::set<dof_id_type> DoF_RCSet;
2618  DoF_RCSet unexpanded_dofs;
2619 
2620  for (DofConstraints::iterator i = _dof_constraints.begin();
2621  i != _dof_constraints.end(); ++i)
2622  {
2623  unexpanded_dofs.insert(i->first);
2624  }
2625 
2626  typedef std::set<const Node *> Node_RCSet;
2627  Node_RCSet unexpanded_nodes;
2628 
2629 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2630  for (NodeConstraints::iterator i = _node_constraints.begin();
2631  i != _node_constraints.end(); ++i)
2632  {
2633  unexpanded_nodes.insert(i->first);
2634  }
2635 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2636 
2637  // We have to keep recursing while the unexpanded set is
2638  // nonempty on *any* processor
2639  bool unexpanded_set_nonempty = !unexpanded_dofs.empty() ||
2640  !unexpanded_nodes.empty();
2641  this->comm().max(unexpanded_set_nonempty);
2642 
2643  while (unexpanded_set_nonempty)
2644  {
2645  // Let's make sure we don't lose sync in this loop.
2646  parallel_object_only();
2647 
2648  // Request sets
2649  DoF_RCSet dof_request_set;
2650  Node_RCSet node_request_set;
2651 
2652  // Request sets to send to each processor
2653  std::vector<std::vector<dof_id_type> >
2654  requested_dof_ids(this->n_processors()),
2655  requested_node_ids(this->n_processors());
2656 
2657  // And the sizes of each
2658  std::vector<dof_id_type>
2659  dof_ids_on_proc(this->n_processors(), 0),
2660  node_ids_on_proc(this->n_processors(), 0);
2661 
2662  // Fill (and thereby sort and uniq!) the main request sets
2663  for (DoF_RCSet::iterator i = unexpanded_dofs.begin();
2664  i != unexpanded_dofs.end(); ++i)
2665  {
2667  for (DofConstraintRow::iterator j = row.begin();
2668  j != row.end(); ++j)
2669  {
2670  const dof_id_type constraining_dof = j->first;
2671 
2672  // If it's non-local and we haven't already got a
2673  // constraint for it, we might need to ask for one
2674  if (((constraining_dof < this->first_dof()) ||
2675  (constraining_dof >= this->end_dof())) &&
2676  !_dof_constraints.count(constraining_dof))
2677  dof_request_set.insert(constraining_dof);
2678  }
2679  }
2680 
2681 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2682  for (Node_RCSet::iterator i = unexpanded_nodes.begin();
2683  i != unexpanded_nodes.end(); ++i)
2684  {
2685  NodeConstraintRow &row = _node_constraints[*i].first;
2686  for (NodeConstraintRow::iterator j = row.begin();
2687  j != row.end(); ++j)
2688  {
2689  const Node* const node = j->first;
2690  libmesh_assert(node);
2691 
2692  // If it's non-local and we haven't already got a
2693  // constraint for it, we might need to ask for one
2694  if ((node->processor_id() != this->processor_id()) &&
2695  !_node_constraints.count(node))
2696  node_request_set.insert(node);
2697  }
2698  }
2699 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2700 
2701  // Clear the unexpanded constraint sets; we're about to expand
2702  // them
2703  unexpanded_dofs.clear();
2704  unexpanded_nodes.clear();
2705 
2706  // Count requests by processor
2707  processor_id_type proc_id = 0;
2708  for (DoF_RCSet::iterator i = dof_request_set.begin();
2709  i != dof_request_set.end(); ++i)
2710  {
2711  while (*i >= _end_df[proc_id])
2712  proc_id++;
2713  dof_ids_on_proc[proc_id]++;
2714  }
2715 
2716  for (Node_RCSet::iterator i = node_request_set.begin();
2717  i != node_request_set.end(); ++i)
2718  {
2719  libmesh_assert(*i);
2720  libmesh_assert_less ((*i)->processor_id(), this->n_processors());
2721  node_ids_on_proc[(*i)->processor_id()]++;
2722  }
2723 
2724  for (processor_id_type p = 0; p != this->n_processors(); ++p)
2725  {
2726  requested_dof_ids[p].reserve(dof_ids_on_proc[p]);
2727  requested_node_ids[p].reserve(node_ids_on_proc[p]);
2728  }
2729 
2730  // Prepare each processor's request set
2731  proc_id = 0;
2732  for (DoF_RCSet::iterator i = dof_request_set.begin();
2733  i != dof_request_set.end(); ++i)
2734  {
2735  while (*i >= _end_df[proc_id])
2736  proc_id++;
2737  requested_dof_ids[proc_id].push_back(*i);
2738  }
2739 
2740  for (Node_RCSet::iterator i = node_request_set.begin();
2741  i != node_request_set.end(); ++i)
2742  {
2743  requested_node_ids[(*i)->processor_id()].push_back((*i)->id());
2744  }
2745 
2746  // Now request constraint rows from other processors
2747  for (processor_id_type p=1; p != this->n_processors(); ++p)
2748  {
2749  // Trade my requests with processor procup and procdown
2750  processor_id_type procup = (this->processor_id() + p) %
2751  this->n_processors();
2752  processor_id_type procdown = (this->n_processors() +
2753  this->processor_id() - p) %
2754  this->n_processors();
2755  std::vector<dof_id_type> dof_request_to_fill,
2756  node_request_to_fill;
2757  this->comm().send_receive(procup, requested_dof_ids[procup],
2758  procdown, dof_request_to_fill);
2759  this->comm().send_receive(procup, requested_node_ids[procup],
2760  procdown, node_request_to_fill);
2761 
2762  // Fill those requests
2763  std::vector<std::vector<dof_id_type> > dof_row_keys(dof_request_to_fill.size()),
2764  node_row_keys(node_request_to_fill.size());
2765 
2766  // FIXME - this could be an unordered set, given a
2767  // hash<pointers> specialization
2768  std::set<const Node*> nodes_requested;
2769 
2770  std::vector<std::vector<Real> > dof_row_vals(dof_request_to_fill.size()),
2771  node_row_vals(node_request_to_fill.size());
2772  std::vector<Number> dof_row_rhss(dof_request_to_fill.size());
2773  std::vector<Point> node_row_rhss(node_request_to_fill.size());
2774  for (std::size_t i=0; i != dof_request_to_fill.size(); ++i)
2775  {
2776  dof_id_type constrained = dof_request_to_fill[i];
2777  if (_dof_constraints.count(constrained))
2778  {
2779  DofConstraintRow &row = _dof_constraints[constrained];
2780  std::size_t row_size = row.size();
2781  dof_row_keys[i].reserve(row_size);
2782  dof_row_vals[i].reserve(row_size);
2783  for (DofConstraintRow::iterator j = row.begin();
2784  j != row.end(); ++j)
2785  {
2786  dof_row_keys[i].push_back(j->first);
2787  dof_row_vals[i].push_back(j->second);
2788 
2789  // We should never have a 0 constraint
2790  // coefficient; that's implicit via sparse
2791  // constraint storage
2792  libmesh_assert(j->second);
2793  }
2794  DofConstraintValueMap::const_iterator rhsit =
2795  _primal_constraint_values.find(constrained);
2796  dof_row_rhss[i] = (rhsit == _primal_constraint_values.end()) ?
2797  0 : rhsit->second;
2798  }
2799  }
2800 
2801 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2802  for (std::size_t i=0; i != node_request_to_fill.size(); ++i)
2803  {
2804  dof_id_type constrained_id = node_request_to_fill[i];
2805  const Node *constrained_node = mesh.node_ptr(constrained_id);
2806  if (_node_constraints.count(constrained_node))
2807  {
2808  const NodeConstraintRow &row = _node_constraints[constrained_node].first;
2809  std::size_t row_size = row.size();
2810  node_row_keys[i].reserve(row_size);
2811  node_row_vals[i].reserve(row_size);
2812  for (NodeConstraintRow::const_iterator j = row.begin();
2813  j != row.end(); ++j)
2814  {
2815  const Node* node = j->first;
2816  node_row_keys[i].push_back(node->id());
2817  node_row_vals[i].push_back(j->second);
2818 
2819  // If we're not sure whether our send
2820  // destination already has this node, let's give
2821  // it a copy.
2822  if (node->processor_id() != procdown)
2823  nodes_requested.insert(node);
2824 
2825  // We can have 0 nodal constraint
2826  // coefficients, where no Lagrange constrant
2827  // exists but non-Lagrange basis constraints
2828  // might.
2829  // libmesh_assert(j->second);
2830  }
2831  node_row_rhss[i] = _node_constraints[constrained_node].second;
2832  }
2833  }
2834 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2835 
2836  // Trade back the results
2837  std::vector<std::vector<dof_id_type> > dof_filled_keys,
2838  node_filled_keys;
2839  std::vector<std::vector<Real> > dof_filled_vals,
2840  node_filled_vals;
2841  std::vector<Number> dof_filled_rhss;
2842  std::vector<Point> node_filled_rhss;
2843  this->comm().send_receive(procdown, dof_row_keys,
2844  procup, dof_filled_keys);
2845  this->comm().send_receive(procdown, dof_row_vals,
2846  procup, dof_filled_vals);
2847  this->comm().send_receive(procdown, dof_row_rhss,
2848  procup, dof_filled_rhss);
2849 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2850  this->comm().send_receive(procdown, node_row_keys,
2851  procup, node_filled_keys);
2852  this->comm().send_receive(procdown, node_row_vals,
2853  procup, node_filled_vals);
2854  this->comm().send_receive(procdown, node_row_rhss,
2855  procup, node_filled_rhss);
2856 
2857  // Constraining nodes might not even exist on our subset of
2858  // a distributed mesh, so let's make them exist.
2860  (procdown, &mesh, nodes_requested.begin(), nodes_requested.end(),
2862 
2863 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2864  libmesh_assert_equal_to (dof_filled_keys.size(), requested_dof_ids[procup].size());
2865  libmesh_assert_equal_to (dof_filled_vals.size(), requested_dof_ids[procup].size());
2866  libmesh_assert_equal_to (dof_filled_rhss.size(), requested_dof_ids[procup].size());
2867  libmesh_assert_equal_to (node_filled_keys.size(), requested_node_ids[procup].size());
2868  libmesh_assert_equal_to (node_filled_vals.size(), requested_node_ids[procup].size());
2869  libmesh_assert_equal_to (node_filled_rhss.size(), requested_node_ids[procup].size());
2870 
2871  // Add any new constraint rows we've found
2872  for (std::size_t i=0; i != requested_dof_ids[procup].size(); ++i)
2873  {
2874  libmesh_assert_equal_to (dof_filled_keys[i].size(), dof_filled_vals[i].size());
2875  // FIXME - what about empty p constraints!?
2876  if (!dof_filled_keys[i].empty())
2877  {
2878  dof_id_type constrained = requested_dof_ids[procup][i];
2879  DofConstraintRow &row = _dof_constraints[constrained];
2880  for (std::size_t j = 0; j != dof_filled_keys[i].size(); ++j)
2881  row[dof_filled_keys[i][j]] = dof_filled_vals[i][j];
2882  if (dof_filled_rhss[i] != Number(0))
2883  _primal_constraint_values[constrained] = dof_filled_rhss[i];
2884  else
2885  _primal_constraint_values.erase(constrained);
2886 
2887  // And prepare to check for more recursive constraints
2888  unexpanded_dofs.insert(constrained);
2889  }
2890  }
2891 
2892 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2893  for (std::size_t i=0; i != requested_node_ids[procup].size(); ++i)
2894  {
2895  libmesh_assert_equal_to (node_filled_keys[i].size(), node_filled_vals[i].size());
2896  if (!node_filled_keys[i].empty())
2897  {
2898  dof_id_type constrained_id = requested_node_ids[procup][i];
2899  const Node* constrained_node = mesh.node_ptr(constrained_id);
2900  NodeConstraintRow &row = _node_constraints[constrained_node].first;
2901  for (std::size_t j = 0; j != node_filled_keys[i].size(); ++j)
2902  {
2903  const Node* key_node =
2904  mesh.node_ptr(node_filled_keys[i][j]);
2905  libmesh_assert(key_node);
2906  row[key_node] = node_filled_vals[i][j];
2907  }
2908  _node_constraints[constrained_node].second = node_filled_rhss[i];
2909 
2910  // And prepare to check for more recursive constraints
2911  unexpanded_nodes.insert(constrained_node);
2912  }
2913  }
2914 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2915  }
2916 
2917  // We have to keep recursing while the unexpanded set is
2918  // nonempty on *any* processor
2919  unexpanded_set_nonempty = !unexpanded_dofs.empty() ||
2920  !unexpanded_nodes.empty();
2921  this->comm().max(unexpanded_set_nonempty);
2922  }
2923 }
void libMesh::DofMap::attach_extra_send_list_function ( void(*)(std::vector< dof_id_type > &, void *)  func,
void *  context = NULL 
)
inline

Attach a function pointer to use as a callback to populate the send_list with extra entries.

Definition at line 298 of file dof_map.h.

References _extra_send_list_context, and _extra_send_list_function.

void libMesh::DofMap::attach_extra_send_list_object ( DofMap::AugmentSendList asl)
inline

Attach an object to populate the send_list with extra entries. This should only add to the send list, but no checking is done to enforce this behavior.

This is an advanced function... use at your own peril!

Definition at line 289 of file dof_map.h.

References _augment_send_list.

290  {
291  _augment_send_list = &asl;
292  }
void libMesh::DofMap::attach_extra_sparsity_function ( void(*)(SparsityPattern::Graph &sparsity, std::vector< dof_id_type > &n_nz, std::vector< dof_id_type > &n_oz, void *)  func,
void *  context = NULL 
)
inline

Attach a function pointer to use as a callback to populate the sparsity pattern with extra entries.

Care must be taken that when adding entries they are sorted into the Rows

Further, you must modify n_nz and n_oz properly!

This is an advanced function... use at your own peril!

Definition at line 275 of file dof_map.h.

References _extra_sparsity_context, and _extra_sparsity_function.

void libMesh::DofMap::attach_extra_sparsity_object ( DofMap::AugmentSparsityPattern asp)
inline

Attach an object to use to populate the sparsity pattern with extra entries.

Care must be taken that when adding entries they are sorted into the Rows

Further, you must modify n_nz and n_oz properly!

This is an advanced function... use at your own peril!

Definition at line 260 of file dof_map.h.

References _augment_sparsity_pattern.

261  {
263  }
void libMesh::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 225 of file dof_map.C.

References _matrices, _n_nz, _n_oz, _sp, libMesh::SparseMatrix< T >::attach_dof_map(), libMesh::ParallelObject::comm(), libMesh::AutoPtr< Tp >::get(), libMesh::libmesh_assert(), libMesh::Parallel::Communicator::max(), libMesh::SparseMatrix< T >::need_full_sparsity_pattern(), need_full_sparsity_pattern, libMesh::SparsityPattern::Build::sparsity_pattern, and libMesh::SparseMatrix< T >::update_sparsity_pattern().

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

226 {
227  parallel_object_only();
228 
229  // We shouldn't be trying to re-attach the same matrices repeatedly
230  libmesh_assert (std::find(_matrices.begin(), _matrices.end(),
231  &matrix) == _matrices.end());
232 
233  _matrices.push_back(&matrix);
234 
235  matrix.attach_dof_map (*this);
236 
237  // If we've already computed sparsity, then it's too late
238  // to wait for "compute_sparsity" to help with sparse matrix
239  // initialization, and we need to handle this matrix individually
240  bool computed_sparsity_already =
241  ((_n_nz && !_n_nz->empty()) ||
242  (_n_oz && !_n_oz->empty()));
243  this->comm().max(computed_sparsity_already);
244  if (computed_sparsity_already &&
245  matrix.need_full_sparsity_pattern())
246  {
247  // We'd better have already computed the full sparsity pattern
248  // if we need it here
251 
252  matrix.update_sparsity_pattern (_sp->sparsity_pattern);
253  }
254 
255  if (matrix.need_full_sparsity_pattern())
257 }
unsigned int libMesh::DofMap::block_size ( ) const
inline
Returns
the block size, if the variables are amenable to block storage. Otherwise 1.

Definition at line 416 of file dof_map.h.

References has_blocked_representation(), and n_variables().

417  {
418 #ifdef LIBMESH_ENABLE_BLOCKED_STORAGE
419  return (this->has_blocked_representation() ? this->n_variables() : 1);
420 #else
421  return 1;
422 #endif
423  }
void libMesh::DofMap::build_constraint_matrix ( DenseMatrix< Number > &  C,
std::vector< dof_id_type > &  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 2134 of file dof_map_constraints.C.

References libMesh::libmesh_assert(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), libMesh::DenseMatrix< T >::resize(), libMesh::DenseMatrix< T >::right_multiply(), libMesh::START_LOG(), and libMesh::STOP_LOG().

2137 {
2138  if (!called_recursively) START_LOG("build_constraint_matrix()", "DofMap");
2139 
2140  // Create a set containing the DOFs we already depend on
2141  typedef std::set<dof_id_type> RCSet;
2142  RCSet dof_set;
2143 
2144  bool we_have_constraints = false;
2145 
2146  // Next insert any other dofs the current dofs might be constrained
2147  // in terms of. Note that in this case we may not be done: Those
2148  // may in turn depend on others. So, we need to repeat this process
2149  // in that case until the system depends only on unconstrained
2150  // degrees of freedom.
2151  for (unsigned int i=0; i<elem_dofs.size(); i++)
2152  if (this->is_constrained_dof(elem_dofs[i]))
2153  {
2154  we_have_constraints = true;
2155 
2156  // If the DOF is constrained
2157  DofConstraints::const_iterator
2158  pos = _dof_constraints.find(elem_dofs[i]);
2159 
2160  libmesh_assert (pos != _dof_constraints.end());
2161 
2162  const DofConstraintRow& constraint_row = pos->second;
2163 
2164 // Constraint rows in p refinement may be empty
2165 // libmesh_assert (!constraint_row.empty());
2166 
2167  for (DofConstraintRow::const_iterator
2168  it=constraint_row.begin(); it != constraint_row.end();
2169  ++it)
2170  dof_set.insert (it->first);
2171  }
2172 
2173  // May be safe to return at this point
2174  // (but remember to stop the perflog)
2175  if (!we_have_constraints)
2176  {
2177  STOP_LOG("build_constraint_matrix()", "DofMap");
2178  return;
2179  }
2180 
2181  for (unsigned int i=0; i != elem_dofs.size(); ++i)
2182  dof_set.erase (elem_dofs[i]);
2183 
2184  // If we added any DOFS then we need to do this recursively.
2185  // It is possible that we just added a DOF that is also
2186  // constrained!
2187  //
2188  // Also, we need to handle the special case of an element having DOFs
2189  // constrained in terms of other, local DOFs
2190  if (!dof_set.empty() || // case 1: constrained in terms of other DOFs
2191  !called_recursively) // case 2: constrained in terms of our own DOFs
2192  {
2193  const unsigned int old_size =
2194  libmesh_cast_int<unsigned int>(elem_dofs.size());
2195 
2196  // Add new dependency dofs to the end of the current dof set
2197  elem_dofs.insert(elem_dofs.end(),
2198  dof_set.begin(), dof_set.end());
2199 
2200  // Now we can build the constraint matrix.
2201  // Note that resize also zeros for a DenseMatrix<Number>.
2202  C.resize (old_size,
2203  libmesh_cast_int<unsigned int>(elem_dofs.size()));
2204 
2205  // Create the C constraint matrix.
2206  for (unsigned int i=0; i != old_size; i++)
2207  if (this->is_constrained_dof(elem_dofs[i]))
2208  {
2209  // If the DOF is constrained
2210  DofConstraints::const_iterator
2211  pos = _dof_constraints.find(elem_dofs[i]);
2212 
2213  libmesh_assert (pos != _dof_constraints.end());
2214 
2215  const DofConstraintRow& constraint_row = pos->second;
2216 
2217 // p refinement creates empty constraint rows
2218 // libmesh_assert (!constraint_row.empty());
2219 
2220  for (DofConstraintRow::const_iterator
2221  it=constraint_row.begin(); it != constraint_row.end();
2222  ++it)
2223  for (unsigned int j=0; j != elem_dofs.size(); j++)
2224  if (elem_dofs[j] == it->first)
2225  C(i,j) = it->second;
2226  }
2227  else
2228  {
2229  C(i,i) = 1.;
2230  }
2231 
2232  // May need to do this recursively. It is possible
2233  // that we just replaced a constrained DOF with another
2234  // constrained DOF.
2235  DenseMatrix<Number> Cnew;
2236 
2237  this->build_constraint_matrix (Cnew, elem_dofs, true);
2238 
2239  if ((C.n() == Cnew.m()) &&
2240  (Cnew.n() == elem_dofs.size())) // If the constraint matrix
2241  C.right_multiply(Cnew); // is constrained...
2242 
2243  libmesh_assert_equal_to (C.n(), elem_dofs.size());
2244  }
2245 
2246  if (!called_recursively) STOP_LOG("build_constraint_matrix()", "DofMap");
2247 }
void libMesh::DofMap::build_constraint_matrix_and_vector ( DenseMatrix< Number > &  C,
DenseVector< Number > &  H,
std::vector< dof_id_type > &  elem_dofs,
int  qoi_index = -1,
const bool  called_recursively = false 
) const
private

Build the constraint matrix C and the forcing vector H 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 and/or Dirichlet conditions.

The forcing vector will depend on which solution's heterogenous constraints are being applied. For the default qoi_index this will be the primal solutoin; for qoi_index >= 0 the corresponding adjoint solution's constraints will be used.

Definition at line 2252 of file dof_map_constraints.C.

References libMesh::libmesh_assert(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::DenseMatrix< T >::right_multiply(), libMesh::START_LOG(), libMesh::STOP_LOG(), and libMesh::DenseMatrix< T >::vector_mult_add().

Referenced by extract_local_vector().

2257 {
2258  if (!called_recursively)
2259  START_LOG("build_constraint_matrix_and_vector()", "DofMap");
2260 
2261  // Create a set containing the DOFs we already depend on
2262  typedef std::set<dof_id_type> RCSet;
2263  RCSet dof_set;
2264 
2265  bool we_have_constraints = false;
2266 
2267  // Next insert any other dofs the current dofs might be constrained
2268  // in terms of. Note that in this case we may not be done: Those
2269  // may in turn depend on others. So, we need to repeat this process
2270  // in that case until the system depends only on unconstrained
2271  // degrees of freedom.
2272  for (unsigned int i=0; i<elem_dofs.size(); i++)
2273  if (this->is_constrained_dof(elem_dofs[i]))
2274  {
2275  we_have_constraints = true;
2276 
2277  // If the DOF is constrained
2278  DofConstraints::const_iterator
2279  pos = _dof_constraints.find(elem_dofs[i]);
2280 
2281  libmesh_assert (pos != _dof_constraints.end());
2282 
2283  const DofConstraintRow& constraint_row = pos->second;
2284 
2285 // Constraint rows in p refinement may be empty
2286 // libmesh_assert (!constraint_row.empty());
2287 
2288  for (DofConstraintRow::const_iterator
2289  it=constraint_row.begin(); it != constraint_row.end();
2290  ++it)
2291  dof_set.insert (it->first);
2292  }
2293 
2294  // May be safe to return at this point
2295  // (but remember to stop the perflog)
2296  if (!we_have_constraints)
2297  {
2298  STOP_LOG("build_constraint_matrix_and_vector()", "DofMap");
2299  return;
2300  }
2301 
2302  for (unsigned int i=0; i != elem_dofs.size(); ++i)
2303  dof_set.erase (elem_dofs[i]);
2304 
2305  // If we added any DOFS then we need to do this recursively.
2306  // It is possible that we just added a DOF that is also
2307  // constrained!
2308  //
2309  // Also, we need to handle the special case of an element having DOFs
2310  // constrained in terms of other, local DOFs
2311  if (!dof_set.empty() || // case 1: constrained in terms of other DOFs
2312  !called_recursively) // case 2: constrained in terms of our own DOFs
2313  {
2314  const DofConstraintValueMap *rhs_values = NULL;
2315  if (qoi_index < 0)
2316  rhs_values = &_primal_constraint_values;
2317  else
2318  {
2319  const AdjointDofConstraintValues::const_iterator
2320  it = _adjoint_constraint_values.find(qoi_index);
2321  if (it != _adjoint_constraint_values.end())
2322  rhs_values = &it->second;
2323  }
2324 
2325  const unsigned int old_size =
2326  libmesh_cast_int<unsigned int>(elem_dofs.size());
2327 
2328  // Add new dependency dofs to the end of the current dof set
2329  elem_dofs.insert(elem_dofs.end(),
2330  dof_set.begin(), dof_set.end());
2331 
2332  // Now we can build the constraint matrix and vector.
2333  // Note that resize also zeros for a DenseMatrix and DenseVector
2334  C.resize (old_size,
2335  libmesh_cast_int<unsigned int>(elem_dofs.size()));
2336  H.resize (old_size);
2337 
2338  // Create the C constraint matrix.
2339  for (unsigned int i=0; i != old_size; i++)
2340  if (this->is_constrained_dof(elem_dofs[i]))
2341  {
2342  // If the DOF is constrained
2343  DofConstraints::const_iterator
2344  pos = _dof_constraints.find(elem_dofs[i]);
2345 
2346  libmesh_assert (pos != _dof_constraints.end());
2347 
2348  const DofConstraintRow& constraint_row = pos->second;
2349 
2350 // p refinement creates empty constraint rows
2351 // libmesh_assert (!constraint_row.empty());
2352 
2353  for (DofConstraintRow::const_iterator
2354  it=constraint_row.begin(); it != constraint_row.end();
2355  ++it)
2356  for (unsigned int j=0; j != elem_dofs.size(); j++)
2357  if (elem_dofs[j] == it->first)
2358  C(i,j) = it->second;
2359 
2360  if (rhs_values)
2361  {
2362  DofConstraintValueMap::const_iterator rhsit =
2363  rhs_values->find(elem_dofs[i]);
2364  if (rhsit != rhs_values->end())
2365  H(i) = rhsit->second;
2366  }
2367  }
2368  else
2369  {
2370  C(i,i) = 1.;
2371  }
2372 
2373  // May need to do this recursively. It is possible
2374  // that we just replaced a constrained DOF with another
2375  // constrained DOF.
2376  DenseMatrix<Number> Cnew;
2377  DenseVector<Number> Hnew;
2378 
2379  this->build_constraint_matrix_and_vector (Cnew, Hnew, elem_dofs,
2380  qoi_index, true);
2381 
2382  if ((C.n() == Cnew.m()) && // If the constraint matrix
2383  (Cnew.n() == elem_dofs.size())) // is constrained...
2384  {
2385  C.right_multiply(Cnew);
2386 
2387  // If x = Cy + h and y = Dz + g
2388  // Then x = (CD)z + (Cg + h)
2389  C.vector_mult_add(H, 1, Hnew);
2390  }
2391 
2392  libmesh_assert_equal_to (C.n(), elem_dofs.size());
2393  }
2394 
2395  if (!called_recursively)
2396  STOP_LOG("build_constraint_matrix_and_vector()", "DofMap");
2397 }
AutoPtr< SparsityPattern::Build > libMesh::DofMap::build_sparsity ( const MeshBase mesh) const
private

Builds a sparsity pattern

Definition at line 54 of file dof_map.C.

References libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::MeshBase::is_prepared(), libMesh::libmesh_assert(), libMesh::SparsityPattern::Build::n_nz, libMesh::SparsityPattern::Build::n_oz, libMesh::out, libMesh::Threads::parallel_reduce(), libMesh::SparsityPattern::Build::parallel_sync(), libMesh::ParallelObject::processor_id(), libMesh::SparsityPattern::Build::sparsity_pattern, libMesh::START_LOG(), and libMesh::STOP_LOG().

Referenced by compute_sparsity().

55 {
56  libmesh_assert (mesh.is_prepared());
57  libmesh_assert (this->n_variables());
58 
59  START_LOG("build_sparsity()", "DofMap");
60 
61  // Compute the sparsity structure of the global matrix. This can be
62  // fed into a PetscMatrix to allocate exacly the number of nonzeros
63  // necessary to store the matrix. This algorithm should be linear
64  // in the (# of elements)*(# nodes per element)
65 
66  // We can be more efficient in the threaded sparsity pattern assembly
67  // if we don't need the exact pattern. For some sparse matrix formats
68  // a good upper bound will suffice.
69 
70  // See if we need to include sparsity pattern entries for coupling
71  // between neighbor dofs
72  bool implicit_neighbor_dofs = this->use_coupled_neighbor_dofs(mesh);
73 
74  // We can compute the sparsity pattern in parallel on multiple
75  // threads. The goal is for each thread to compute the full sparsity
76  // pattern for a subset of elements. These sparsity patterns can
77  // be efficiently merged in the SparsityPattern::Build::join()
78  // method, especially if there is not too much overlap between them.
79  // Even better, if the full sparsity pattern is not needed then
80  // the number of nonzeros per row can be estimated from the
81  // sparsity patterns created on each thread.
82  AutoPtr<SparsityPattern::Build> sp
83  (new SparsityPattern::Build (mesh,
84  *this,
85  this->_dof_coupling,
86  implicit_neighbor_dofs,
88 
89  Threads::parallel_reduce (ConstElemRange (mesh.active_local_elements_begin(),
90  mesh.active_local_elements_end()), *sp);
91 
92  sp->parallel_sync();
93 
94 #ifndef NDEBUG
95  // Avoid declaring these variables unless asserts are enabled.
96  const processor_id_type proc_id = mesh.processor_id();
97  const dof_id_type n_dofs_on_proc = this->n_dofs_on_processor(proc_id);
98 #endif
99  libmesh_assert_equal_to (sp->sparsity_pattern.size(), n_dofs_on_proc);
100 
101  STOP_LOG("build_sparsity()", "DofMap");
102 
103  // Check to see if we have any extra stuff to add to the sparsity_pattern
105  {
107  {
108  libmesh_here();
109  libMesh::out << "WARNING: You have specified both an extra sparsity function and object.\n"
110  << " Are you sure this is what you meant to do??"
111  << std::endl;
112  }
113 
115  (sp->sparsity_pattern, sp->n_nz,
116  sp->n_oz, _extra_sparsity_context);
117  }
118 
121  (sp->sparsity_pattern, sp->n_nz, sp->n_oz);
122 
123  return sp;
124 }
void libMesh::DofMap::clear ( )

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

Definition at line 808 of file dof_map.C.

References _adjoint_constraint_values, _dof_constraints, _end_df, _end_old_df, _first_df, _first_old_df, _matrices, _n_dfs, _n_old_dfs, _primal_constraint_values, _send_list, _variable_groups, _variables, clear_sparsity(), and need_full_sparsity_pattern.

Referenced by DofMap().

809 {
810  // we don't want to clear
811  // the coupling matrix!
812  // It should not change...
813  //_dof_coupling->clear();
814 
815  _variables.clear();
816  _variable_groups.clear();
817  _first_df.clear();
818  _end_df.clear();
819  _send_list.clear();
820  this->clear_sparsity();
822 
823 #ifdef LIBMESH_ENABLE_AMR
824 
825  _dof_constraints.clear();
828  _n_old_dfs = 0;
829  _first_old_df.clear();
830  _end_old_df.clear();
831 
832 #endif
833 
834  _matrices.clear();
835 
836  _n_dfs = 0;
837 }
void libMesh::DofMap::clear_sparsity ( )

Clears the sparsity pattern

Definition at line 1531 of file dof_map.C.

References _n_nz, _n_oz, _sp, libMesh::AutoPtr< Tp >::get(), libMesh::libmesh_assert(), libMesh::SparsityPattern::Build::n_nz, libMesh::SparsityPattern::Build::n_oz, need_full_sparsity_pattern, and libMesh::AutoPtr< Tp >::reset().

Referenced by clear(), libMesh::ImplicitSystem::reinit(), and libMesh::EigenSystem::reinit().

1532 {
1534  {
1535  libmesh_assert(_sp.get());
1536  libmesh_assert(!_n_nz || _n_nz == &_sp->n_nz);
1537  libmesh_assert(!_n_oz || _n_oz == &_sp->n_oz);
1538  _sp.reset();
1539  }
1540  else
1541  {
1542  libmesh_assert(!_sp.get());
1543  delete _n_nz;
1544  delete _n_oz;
1545  }
1546  _n_nz = NULL;
1547  _n_oz = NULL;
1548 }
const Parallel::Communicator& libMesh::ParallelObject::comm ( ) const
inlineinherited
Returns
a reference to the Parallel::Communicator object used by this mesh.

Definition at line 86 of file parallel_object.h.

References libMesh::ParallelObject::_communicator.

Referenced by libMesh::__libmesh_petsc_diff_solver_jacobian(), libMesh::__libmesh_petsc_diff_solver_monitor(), libMesh::__libmesh_petsc_diff_solver_residual(), libMesh::__libmesh_petsc_snes_jacobian(), libMesh::__libmesh_petsc_snes_residual(), libMesh::MeshRefinement::_coarsen_elements(), libMesh::ExactSolution::_compute_error(), libMesh::MetisPartitioner::_do_partition(), libMesh::ParmetisPartitioner::_do_repartition(), libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_get_diagonal(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_get_diagonal(), libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_mult(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult_add(), libMesh::EquationSystems::_read_impl(), libMesh::MeshRefinement::_refine_elements(), libMesh::ParallelMesh::add_elem(), libMesh::ImplicitSystem::add_matrix(), libMesh::ParallelMesh::add_node(), libMesh::System::add_vector(), libMesh::UnstructuredMesh::all_second_order(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::FEMSystem::assemble_qoi(), libMesh::MeshCommunication::assign_global_indices(), libMesh::ParmetisPartitioner::assign_partitioning(), attach_matrix(), libMesh::MeshTools::bounding_box(), libMesh::System::calculate_norm(), libMesh::MeshRefinement::coarsen_elements(), libMesh::Nemesis_IO_Helper::compute_num_global_elem_blocks(), libMesh::Nemesis_IO_Helper::compute_num_global_nodesets(), libMesh::Nemesis_IO_Helper::compute_num_global_sidesets(), libMesh::Problem_Interface::computeF(), libMesh::Problem_Interface::computeJacobian(), libMesh::Problem_Interface::computePreconditioner(), libMesh::MeshTools::correct_node_proc_ids(), libMesh::MeshCommunication::delete_remote_elements(), distribute_dofs(), DMlibMeshFunction(), DMlibMeshJacobian(), DMLibMeshSetSystem(), DMVariableBounds_libMesh(), libMesh::MeshRefinement::eliminate_unrefined_patches(), libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::MeshRefinement::flag_elements_by_elem_fraction(), libMesh::MeshRefinement::flag_elements_by_error_fraction(), libMesh::MeshRefinement::flag_elements_by_nelem_target(), libMesh::for(), libMesh::CondensedEigenSystem::get_eigenpair(), libMesh::ImplicitSystem::get_linear_solver(), libMesh::LocationMap< T >::init(), libMesh::PetscDiffSolver::init(), libMesh::TimeSolver::init(), libMesh::SystemSubsetBySubdomain::init(), libMesh::EigenSystem::init_data(), libMesh::EigenSystem::init_matrices(), libMesh::ParmetisPartitioner::initialize(), libMesh::MeshTools::libmesh_assert_valid_dof_ids(), libMesh::ParallelMesh::libmesh_assert_valid_parallel_flags(), libMesh::MeshTools::libmesh_assert_valid_procids< Elem >(), libMesh::MeshTools::libmesh_assert_valid_procids< Node >(), libMesh::MeshTools::libmesh_assert_valid_refinement_flags(), libMesh::MeshRefinement::limit_level_mismatch_at_edge(), libMesh::MeshRefinement::limit_level_mismatch_at_node(), libMesh::MeshRefinement::make_coarsening_compatible(), libMesh::MeshCommunication::make_elems_parallel_consistent(), libMesh::MeshRefinement::make_flags_parallel_consistent(), libMesh::MeshCommunication::make_node_ids_parallel_consistent(), libMesh::MeshCommunication::make_node_proc_ids_parallel_consistent(), libMesh::MeshCommunication::make_nodes_parallel_consistent(), libMesh::MeshRefinement::make_refinement_compatible(), libMesh::FEMSystem::mesh_position_set(), libMesh::MeshSerializer::MeshSerializer(), libMesh::ParallelMesh::n_active_elem(), libMesh::MeshTools::n_active_levels(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::BoundaryInfo::n_edge_conds(), libMesh::CondensedEigenSystem::n_global_non_condensed_dofs(), libMesh::MeshTools::n_levels(), libMesh::BoundaryInfo::n_nodeset_conds(), libMesh::MeshTools::n_p_levels(), libMesh::ParallelMesh::parallel_max_elem_id(), libMesh::ParallelMesh::parallel_max_node_id(), libMesh::ParallelMesh::parallel_n_elem(), libMesh::ParallelMesh::parallel_n_nodes(), libMesh::Partitioner::partition(), libMesh::Partitioner::partition_unpartitioned_elements(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::MeshBase::prepare_for_use(), libMesh::System::project_vector(), libMesh::Nemesis_IO::read(), libMesh::XdrIO::read(), libMesh::System::read_header(), libMesh::System::read_legacy_data(), libMesh::System::read_SCALAR_dofs(), libMesh::XdrIO::read_serialized_bc_names(), libMesh::XdrIO::read_serialized_bcs(), libMesh::System::read_serialized_blocked_dof_objects(), libMesh::XdrIO::read_serialized_connectivity(), libMesh::XdrIO::read_serialized_nodes(), libMesh::XdrIO::read_serialized_nodesets(), libMesh::XdrIO::read_serialized_subdomain_names(), libMesh::System::read_serialized_vector(), libMesh::MeshBase::recalculate_n_partitions(), libMesh::MeshRefinement::refine_and_coarsen_elements(), libMesh::MeshRefinement::refine_elements(), libMesh::Partitioner::set_node_processor_ids(), set_nonlocal_dof_objects(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::MeshBase::subdomain_ids(), libMesh::BoundaryInfo::sync(), libMesh::Parallel::sync_element_data_by_parent_id(), libMesh::MeshRefinement::test_level_one(), libMesh::MeshRefinement::test_unflagged(), libMesh::MeshTools::total_weight(), libMesh::CheckpointIO::write(), libMesh::XdrIO::write(), libMesh::UnstructuredMesh::write(), libMesh::LegacyXdrIO::write_mesh(), libMesh::System::write_SCALAR_dofs(), libMesh::XdrIO::write_serialized_bcs(), libMesh::System::write_serialized_blocked_dof_objects(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), libMesh::XdrIO::write_serialized_nodesets(), and libMesh::DivaIO::write_stream().

87  { return _communicator; }
void libMesh::DofMap::compute_sparsity ( const MeshBase mesh)

Computes the sparsity pattern for the matrices corresponding to proc_id and sends that data to Linear Algebra packages for preallocation of sparse matrices.

Definition at line 1494 of file dof_map.C.

References _matrices, _n_nz, _n_oz, _sp, build_sparsity(), end, libMesh::SparsityPattern::Build::n_nz, libMesh::SparsityPattern::Build::n_oz, need_full_sparsity_pattern, libMesh::AutoPtr< Tp >::reset(), and libMesh::SparsityPattern::Build::sparsity_pattern.

Referenced by libMesh::EigenSystem::init_matrices(), libMesh::ImplicitSystem::init_matrices(), libMesh::ImplicitSystem::reinit(), and libMesh::EigenSystem::reinit().

1495 {
1496  _sp = this->build_sparsity(mesh);
1497 
1498  // It is possible that some \p SparseMatrix implementations want to
1499  // see it. Let them see it before we throw it away.
1500  std::vector<SparseMatrix<Number>* >::const_iterator
1501  pos = _matrices.begin(),
1502  end = _matrices.end();
1503 
1504  // If we need the full sparsity pattern, then we share a view of its
1505  // arrays, and we pass it in to the matrices.
1507  {
1508  _n_nz = &_sp->n_nz;
1509  _n_oz = &_sp->n_oz;
1510 
1511  for (; pos != end; ++pos)
1512  (*pos)->update_sparsity_pattern (_sp->sparsity_pattern);
1513  }
1514  // If we don't need the full sparsity pattern anymore, steal the
1515  // arrays we do need and free the rest of the memory
1516  else
1517  {
1518  if (!_n_nz)
1519  _n_nz = new std::vector<dof_id_type>();
1520  _n_nz->swap(_sp->n_nz);
1521  if (!_n_oz)
1522  _n_oz = new std::vector<dof_id_type>();
1523  _n_oz->swap(_sp->n_oz);
1524 
1525  _sp.reset();
1526  }
1527 }
void libMesh::DofMap::constrain_element_dyad_matrix ( DenseVector< Number > &  v,
DenseVector< Number > &  w,
std::vector< dof_id_type > &  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 1770 of file dof_map_constraints.C.

References libMesh::libmesh_assert(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), libMesh::DenseVector< T >::size(), libMesh::START_LOG(), libMesh::STOP_LOG(), and libMesh::DenseMatrix< T >::vector_mult_transpose().

1774 {
1775  libmesh_assert_equal_to (v.size(), row_dofs.size());
1776  libmesh_assert_equal_to (w.size(), row_dofs.size());
1777 
1778  // check for easy return
1779  if (this->_dof_constraints.empty())
1780  return;
1781 
1782  // The constrained RHS is built up as R^T F.
1784 
1785  this->build_constraint_matrix (R, row_dofs);
1786 
1787  START_LOG("cnstrn_elem_dyad_mat()", "DofMap");
1788 
1789  // It is possible that the vector is not constrained at all.
1790  if ((R.m() == v.size()) &&
1791  (R.n() == row_dofs.size())) // if the RHS is constrained
1792  {
1793  // Compute the matrix-vector products
1794  DenseVector<Number> old_v(v);
1795  DenseVector<Number> old_w(w);
1796 
1797  // compute matrix/vector product
1798  R.vector_mult_transpose(v, old_v);
1799  R.vector_mult_transpose(w, old_w);
1800 
1801  libmesh_assert_equal_to (row_dofs.size(), v.size());
1802  libmesh_assert_equal_to (row_dofs.size(), w.size());
1803 
1804  /* Constrain only v, not w. */
1805 
1806  for (unsigned int i=0; i<row_dofs.size(); i++)
1807  if (this->is_constrained_dof(row_dofs[i]))
1808  {
1809  // If the DOF is constrained
1810  libmesh_assert (_dof_constraints.find(row_dofs[i]) != _dof_constraints.end());
1811 
1812  v(i) = 0;
1813  }
1814  } // end if the RHS is constrained.
1815 
1816  STOP_LOG("cnstrn_elem_dyad_mat()", "DofMap");
1817 }
void libMesh::DofMap::constrain_element_matrix ( DenseMatrix< Number > &  matrix,
std::vector< dof_id_type > &  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 1304 of file dof_map_constraints.C.

References libMesh::DenseMatrix< T >::left_multiply_transpose(), libMesh::libmesh_assert(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), libMesh::DenseMatrix< T >::right_multiply(), libMesh::START_LOG(), and libMesh::STOP_LOG().

1307 {
1308  libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
1309  libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
1310 
1311  // check for easy return
1312  if (this->_dof_constraints.empty())
1313  return;
1314 
1315  // The constrained matrix is built up as C^T K C.
1317 
1318 
1319  this->build_constraint_matrix (C, elem_dofs);
1320 
1321  START_LOG("constrain_elem_matrix()", "DofMap");
1322 
1323  // It is possible that the matrix is not constrained at all.
1324  if ((C.m() == matrix.m()) &&
1325  (C.n() == elem_dofs.size())) // It the matrix is constrained
1326  {
1327  // Compute the matrix-matrix-matrix product C^T K C
1328  matrix.left_multiply_transpose (C);
1329  matrix.right_multiply (C);
1330 
1331 
1332  libmesh_assert_equal_to (matrix.m(), matrix.n());
1333  libmesh_assert_equal_to (matrix.m(), elem_dofs.size());
1334  libmesh_assert_equal_to (matrix.n(), elem_dofs.size());
1335 
1336 
1337  for (unsigned int i=0; i<elem_dofs.size(); i++)
1338  // If the DOF is constrained
1339  if (this->is_constrained_dof(elem_dofs[i]))
1340  {
1341  for (unsigned int j=0; j<matrix.n(); j++)
1342  matrix(i,j) = 0.;
1343 
1344  matrix(i,i) = 1.;
1345 
1346  if (asymmetric_constraint_rows)
1347  {
1348  DofConstraints::const_iterator
1349  pos = _dof_constraints.find(elem_dofs[i]);
1350 
1351  libmesh_assert (pos != _dof_constraints.end());
1352 
1353  const DofConstraintRow& constraint_row = pos->second;
1354 
1355  // This is an overzealous assertion in the presence of
1356  // heterogenous constraints: we now can constrain "u_i = c"
1357  // with no other u_j terms involved.
1358  //
1359  // libmesh_assert (!constraint_row.empty());
1360 
1361  for (DofConstraintRow::const_iterator
1362  it=constraint_row.begin(); it != constraint_row.end();
1363  ++it)
1364  for (unsigned int j=0; j<elem_dofs.size(); j++)
1365  if (elem_dofs[j] == it->first)
1366  matrix(i,j) = -it->second;
1367  }
1368  }
1369  } // end if is constrained...
1370 
1371  STOP_LOG("constrain_elem_matrix()", "DofMap");
1372 }
void libMesh::DofMap::constrain_element_matrix ( DenseMatrix< Number > &  matrix,
std::vector< dof_id_type > &  row_dofs,
std::vector< dof_id_type > &  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 1645 of file dof_map_constraints.C.

References libMesh::DenseMatrix< T >::left_multiply_transpose(), libMesh::libmesh_assert(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), libMesh::DenseMatrix< T >::right_multiply(), libMesh::START_LOG(), and libMesh::STOP_LOG().

1649 {
1650  libmesh_assert_equal_to (row_dofs.size(), matrix.m());
1651  libmesh_assert_equal_to (col_dofs.size(), matrix.n());
1652 
1653  // check for easy return
1654  if (this->_dof_constraints.empty())
1655  return;
1656 
1657  // The constrained matrix is built up as R^T K C.
1660 
1661  // Safeguard against the user passing us the same
1662  // object for row_dofs and col_dofs. If that is done
1663  // the calls to build_matrix would fail
1664  std::vector<dof_id_type> orig_row_dofs(row_dofs);
1665  std::vector<dof_id_type> orig_col_dofs(col_dofs);
1666 
1667  this->build_constraint_matrix (R, orig_row_dofs);
1668  this->build_constraint_matrix (C, orig_col_dofs);
1669 
1670  START_LOG("constrain_elem_matrix()", "DofMap");
1671 
1672  row_dofs = orig_row_dofs;
1673  col_dofs = orig_col_dofs;
1674 
1675 
1676  // It is possible that the matrix is not constrained at all.
1677  if ((R.m() == matrix.m()) &&
1678  (R.n() == row_dofs.size()) &&
1679  (C.m() == matrix.n()) &&
1680  (C.n() == col_dofs.size())) // If the matrix is constrained
1681  {
1682  // K_constrained = R^T K C
1683  matrix.left_multiply_transpose (R);
1684  matrix.right_multiply (C);
1685 
1686 
1687  libmesh_assert_equal_to (matrix.m(), row_dofs.size());
1688  libmesh_assert_equal_to (matrix.n(), col_dofs.size());
1689 
1690 
1691  for (unsigned int i=0; i<row_dofs.size(); i++)
1692  if (this->is_constrained_dof(row_dofs[i]))
1693  {
1694  for (unsigned int j=0; j<matrix.n(); j++)
1695  {
1696  if(row_dofs[i] != col_dofs[j])
1697  matrix(i,j) = 0.;
1698  else // If the DOF is constrained
1699  matrix(i,j) = 1.;
1700  }
1701 
1702  if (asymmetric_constraint_rows)
1703  {
1704  DofConstraints::const_iterator
1705  pos = _dof_constraints.find(row_dofs[i]);
1706 
1707  libmesh_assert (pos != _dof_constraints.end());
1708 
1709  const DofConstraintRow& constraint_row = pos->second;
1710 
1711  libmesh_assert (!constraint_row.empty());
1712 
1713  for (DofConstraintRow::const_iterator
1714  it=constraint_row.begin(); it != constraint_row.end();
1715  ++it)
1716  for (unsigned int j=0; j<col_dofs.size(); j++)
1717  if (col_dofs[j] == it->first)
1718  matrix(i,j) = -it->second;
1719  }
1720  }
1721  } // end if is constrained...
1722 
1723  STOP_LOG("constrain_elem_matrix()", "DofMap");
1724 }
void libMesh::DofMap::constrain_element_matrix_and_vector ( DenseMatrix< Number > &  matrix,
DenseVector< Number > &  rhs,
std::vector< dof_id_type > &  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 1376 of file dof_map_constraints.C.

References libMesh::DenseMatrix< T >::left_multiply_transpose(), libMesh::libmesh_assert(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), libMesh::DenseMatrix< T >::right_multiply(), libMesh::DenseVector< T >::size(), libMesh::START_LOG(), libMesh::STOP_LOG(), and libMesh::DenseMatrix< T >::vector_mult_transpose().

1380 {
1381  libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
1382  libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
1383  libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
1384 
1385  // check for easy return
1386  if (this->_dof_constraints.empty())
1387  return;
1388 
1389  // The constrained matrix is built up as C^T K C.
1390  // The constrained RHS is built up as C^T F
1392 
1393  this->build_constraint_matrix (C, elem_dofs);
1394 
1395  START_LOG("cnstrn_elem_mat_vec()", "DofMap");
1396 
1397  // It is possible that the matrix is not constrained at all.
1398  if ((C.m() == matrix.m()) &&
1399  (C.n() == elem_dofs.size())) // It the matrix is constrained
1400  {
1401  // Compute the matrix-matrix-matrix product C^T K C
1402  matrix.left_multiply_transpose (C);
1403  matrix.right_multiply (C);
1404 
1405 
1406  libmesh_assert_equal_to (matrix.m(), matrix.n());
1407  libmesh_assert_equal_to (matrix.m(), elem_dofs.size());
1408  libmesh_assert_equal_to (matrix.n(), elem_dofs.size());
1409 
1410 
1411  for (unsigned int i=0; i<elem_dofs.size(); i++)
1412  if (this->is_constrained_dof(elem_dofs[i]))
1413  {
1414  for (unsigned int j=0; j<matrix.n(); j++)
1415  matrix(i,j) = 0.;
1416 
1417  // If the DOF is constrained
1418  matrix(i,i) = 1.;
1419 
1420  // This will put a nonsymmetric entry in the constraint
1421  // row to ensure that the linear system produces the
1422  // correct value for the constrained DOF.
1423  if (asymmetric_constraint_rows)
1424  {
1425  DofConstraints::const_iterator
1426  pos = _dof_constraints.find(elem_dofs[i]);
1427 
1428  libmesh_assert (pos != _dof_constraints.end());
1429 
1430  const DofConstraintRow& constraint_row = pos->second;
1431 
1432 // p refinement creates empty constraint rows
1433 // libmesh_assert (!constraint_row.empty());
1434 
1435  for (DofConstraintRow::const_iterator
1436  it=constraint_row.begin(); it != constraint_row.end();
1437  ++it)
1438  for (unsigned int j=0; j<elem_dofs.size(); j++)
1439  if (elem_dofs[j] == it->first)
1440  matrix(i,j) = -it->second;
1441  }
1442  }
1443 
1444 
1445  // Compute the matrix-vector product C^T F
1446  DenseVector<Number> old_rhs(rhs);
1447 
1448  // compute matrix/vector product
1449  C.vector_mult_transpose(rhs, old_rhs);
1450  } // end if is constrained...
1451 
1452  STOP_LOG("cnstrn_elem_mat_vec()", "DofMap");
1453 }
void libMesh::DofMap::constrain_element_vector ( DenseVector< Number > &  rhs,
std::vector< dof_id_type > &  dofs,
bool  asymmetric_constraint_rows = true 
) const
inline

Constrains the element vector.

Definition at line 1728 of file dof_map_constraints.C.

References libMesh::libmesh_assert(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), libMesh::DenseVector< T >::size(), libMesh::START_LOG(), libMesh::STOP_LOG(), and libMesh::DenseMatrix< T >::vector_mult_transpose().

1731 {
1732  libmesh_assert_equal_to (rhs.size(), row_dofs.size());
1733 
1734  // check for easy return
1735  if (this->_dof_constraints.empty())
1736  return;
1737 
1738  // The constrained RHS is built up as R^T F.
1740 
1741  this->build_constraint_matrix (R, row_dofs);
1742 
1743  START_LOG("constrain_elem_vector()", "DofMap");
1744 
1745  // It is possible that the vector is not constrained at all.
1746  if ((R.m() == rhs.size()) &&
1747  (R.n() == row_dofs.size())) // if the RHS is constrained
1748  {
1749  // Compute the matrix-vector product
1750  DenseVector<Number> old_rhs(rhs);
1751  R.vector_mult_transpose(rhs, old_rhs);
1752 
1753  libmesh_assert_equal_to (row_dofs.size(), rhs.size());
1754 
1755  for (unsigned int i=0; i<row_dofs.size(); i++)
1756  if (this->is_constrained_dof(row_dofs[i]))
1757  {
1758  // If the DOF is constrained
1759  libmesh_assert (_dof_constraints.find(row_dofs[i]) != _dof_constraints.end());
1760 
1761  rhs(i) = 0;
1762  }
1763  } // end if the RHS is constrained.
1764 
1765  STOP_LOG("constrain_elem_vector()", "DofMap");
1766 }
void libMesh::DofMap::constrain_nothing ( std::vector< dof_id_type > &  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 1821 of file dof_map_constraints.C.

1822 {
1823  // check for easy return
1824  if (this->_dof_constraints.empty())
1825  return;
1826 
1827  // All the work is done by \p build_constraint_matrix. We just need
1828  // a dummy matrix.
1830  this->build_constraint_matrix (R, dofs);
1831 }
void libMesh::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 3490 of file dof_map_constraints.C.

References libMesh::Elem::dim(), libMesh::DofObject::dof_number(), libMesh::Elem::get_node(), libMesh::Elem::is_node_on_side(), libMesh::Elem::is_vertex(), libMesh::libmesh_assert_greater(), libMesh::DofObject::n_comp(), n_nodes, libMesh::Elem::n_nodes(), libMesh::Elem::n_sides(), libMesh::FEType::order, libMesh::Elem::p_level(), libMesh::Threads::spin_mtx, and libMesh::Elem::type().

Referenced by libMesh::FEGenericBase< T >::compute_periodic_constraints(), and libMesh::FEGenericBase< T >::compute_proj_constraints().

3494 {
3495  // We're constraining dofs on elem which correspond to p refinement
3496  // levels above p - this only makes sense if elem's p refinement
3497  // level is above p.
3498  libmesh_assert_greater (elem->p_level(), p);
3499  libmesh_assert_less (s, elem->n_sides());
3500 
3501  const unsigned int sys_num = this->sys_number();
3502  const unsigned int dim = elem->dim();
3503  ElemType type = elem->type();
3504  FEType low_p_fe_type = this->variable_type(var);
3505  FEType high_p_fe_type = this->variable_type(var);
3506  low_p_fe_type.order = static_cast<Order>(low_p_fe_type.order + p);
3507  high_p_fe_type.order = static_cast<Order>(high_p_fe_type.order +
3508  elem->p_level());
3509 
3510  const unsigned int n_nodes = elem->n_nodes();
3511  for (unsigned int n = 0; n != n_nodes; ++n)
3512  if (elem->is_node_on_side(n, s))
3513  {
3514  const Node * const node = elem->get_node(n);
3515  const unsigned int low_nc =
3516  FEInterface::n_dofs_at_node (dim, low_p_fe_type, type, n);
3517  const unsigned int high_nc =
3518  FEInterface::n_dofs_at_node (dim, high_p_fe_type, type, n);
3519 
3520  // since we may be running this method concurretly
3521  // on multiple threads we need to acquire a lock
3522  // before modifying the _dof_constraints object.
3523  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
3524 
3525  if (elem->is_vertex(n))
3526  {
3527  // Add "this is zero" constraint rows for high p vertex
3528  // dofs
3529  for (unsigned int i = low_nc; i != high_nc; ++i)
3530  {
3531  _dof_constraints[node->dof_number(sys_num,var,i)].clear();
3532  _primal_constraint_values.erase(node->dof_number(sys_num,var,i));
3533  }
3534  }
3535  else
3536  {
3537  const unsigned int total_dofs = node->n_comp(sys_num, var);
3538  libmesh_assert_greater_equal (total_dofs, high_nc);
3539  // Add "this is zero" constraint rows for high p
3540  // non-vertex dofs, which are numbered in reverse
3541  for (unsigned int j = low_nc; j != high_nc; ++j)
3542  {
3543  const unsigned int i = total_dofs - j - 1;
3544  _dof_constraints[node->dof_number(sys_num,var,i)].clear();
3545  _primal_constraint_values.erase(node->dof_number(sys_num,var,i));
3546  }
3547  }
3548  }
3549 }
DofConstraints::const_iterator libMesh::DofMap::constraint_rows_begin ( ) const
inline

Returns an iterator pointing to the first DoF constraint row

Definition at line 646 of file dof_map.h.

References _dof_constraints.

647  { return _dof_constraints.begin(); }
DofConstraints::const_iterator libMesh::DofMap::constraint_rows_end ( ) const
inline

Returns an iterator pointing just past the last DoF constraint row

Definition at line 652 of file dof_map.h.

References _dof_constraints.

653  { return _dof_constraints.end(); }
void libMesh::DofMap::create_dof_constraints ( const MeshBase mesh,
Real  time = 0 
)

Rebuilds the raw degree of freedom and DofObject constraints. A time is specified for use in building time-dependent Dirichlet constraints.

Definition at line 941 of file dof_map_constraints.C.

References libMesh::comm, libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), libMesh::StoredRange< iterator_type, object_type >::empty(), libMesh::MeshBase::is_prepared(), libMesh::MeshBase::is_serial(), libMesh::libmesh_assert(), libMesh::MeshBase::local_elements_begin(), libMesh::MeshBase::local_elements_end(), libMesh::MeshBase::mesh_dimension(), libMesh::Threads::parallel_for(), libMesh::StoredRange< iterator_type, object_type >::reset(), libMesh::START_LOG(), libMesh::STOP_LOG(), libMesh::MeshBase::sub_point_locator(), and libMesh::Parallel::verify().

Referenced by libMesh::EquationSystems::allgather(), and libMesh::EquationSystems::reinit().

942 {
943  parallel_object_only();
944 
945  START_LOG("create_dof_constraints()", "DofMap");
946 
947  libmesh_assert (mesh.is_prepared());
948 
949  const unsigned int dim = mesh.mesh_dimension();
950 
951  // We might get constraint equations from AMR hanging nodes in 2D/3D
952  // or from boundary conditions in any dimension
953  const bool possible_local_constraints = false
954 #ifdef LIBMESH_ENABLE_AMR
955  || dim > 1
956 #endif
957 #ifdef LIBMESH_ENABLE_PERIODIC
958  || !_periodic_boundaries->empty()
959 #endif
960 #ifdef LIBMESH_ENABLE_DIRICHLET
961  || !_dirichlet_boundaries->empty()
962 #endif
963  ;
964 
965  // Even if we don't have constraints, another processor might.
966  bool possible_global_constraints = possible_local_constraints;
967 #if defined(LIBMESH_ENABLE_PERIODIC) || defined(LIBMESH_ENABLE_DIRICHLET) || defined(LIBMESH_ENABLE_AMR)
968  libmesh_assert(this->comm().verify(mesh.is_serial()));
969 
970  if (!mesh.is_serial())
971  this->comm().max(possible_global_constraints);
972 #endif
973 
974  if (!possible_global_constraints)
975  {
976  // Clear out any old constraints; maybe the user just deleted
977  // their last remaining dirichlet/periodic/user constraint?
978 #ifdef LIBMESH_ENABLE_CONSTRAINTS
979  _dof_constraints.clear();
982 #endif
983 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
984  _node_constraints.clear();
985 #endif
986 
987  // make sure we stop logging though
988  STOP_LOG("create_dof_constraints()", "DofMap");
989  return;
990  }
991 
992  // Here we build the hanging node constraints. This is done
993  // by enforcing the condition u_a = u_b along hanging sides.
994  // u_a = u_b is collocated at the nodes of side a, which gives
995  // one row of the constraint matrix.
996 
997  // define the range of elements of interest
998  ConstElemRange range;
999  if (possible_local_constraints)
1000  {
1001  // With SerialMesh or a serial ParallelMesh, every processor
1002  // computes every constraint
1004  elem_begin = mesh.elements_begin(),
1005  elem_end = mesh.elements_end();
1006 
1007  // With a parallel ParallelMesh, processors compute only
1008  // their local constraints
1009  if (!mesh.is_serial())
1010  {
1011  elem_begin = mesh.local_elements_begin();
1012  elem_end = mesh.local_elements_end();
1013  }
1014 
1015  // set the range to contain the specified elements
1016  range.reset (elem_begin, elem_end);
1017  }
1018  else
1019  range.reset (mesh.elements_end(), mesh.elements_end());
1020 
1021  // compute_periodic_constraints requires a point_locator() from our
1022  // Mesh, that point_locator() construction is threaded. Rather than
1023  // nest threads within threads we'll make sure it's preconstructed.
1024 #ifdef LIBMESH_ENABLE_PERIODIC
1025  if (!_periodic_boundaries->empty() && !range.empty())
1026  mesh.sub_point_locator();
1027 #endif
1028 
1029 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1030  // recalculate node constraints from scratch
1031  _node_constraints.clear();
1032 
1033  Threads::parallel_for (range,
1034  ComputeNodeConstraints (_node_constraints,
1035  *this,
1036 #ifdef LIBMESH_ENABLE_PERIODIC
1038 #endif
1039  mesh));
1040 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
1041 
1042 
1043  // recalculate dof constraints from scratch
1044  _dof_constraints.clear();
1045  _primal_constraint_values.clear();
1047 
1048  // Look at all the variables in the system. Reset the element
1049  // range at each iteration -- there is no need to reconstruct it.
1050  for (unsigned int variable_number=0; variable_number<this->n_variables();
1051  ++variable_number, range.reset())
1052  Threads::parallel_for (range,
1053  ComputeConstraints (_dof_constraints,
1054  *this,
1055 #ifdef LIBMESH_ENABLE_PERIODIC
1057 #endif
1058  mesh,
1059  variable_number));
1060 
1061 #ifdef LIBMESH_ENABLE_DIRICHLET
1062  for (DirichletBoundaries::iterator
1063  i = _dirichlet_boundaries->begin();
1064  i != _dirichlet_boundaries->end(); ++i, range.reset())
1065  {
1067  (range,
1068  ConstrainDirichlet(*this, mesh, time, **i,
1069  AddPrimalConstraint(*this))
1070  );
1071  }
1072 
1073  for (unsigned int qoi_index = 0;
1074  qoi_index != _adjoint_dirichlet_boundaries.size();
1075  ++qoi_index)
1076  {
1077  for (DirichletBoundaries::iterator
1078  i = _adjoint_dirichlet_boundaries[qoi_index]->begin();
1079  i != _adjoint_dirichlet_boundaries[qoi_index]->end();
1080  ++i, range.reset())
1081  {
1083  (range,
1084  ConstrainDirichlet(*this, mesh, time, **i,
1085  AddAdjointConstraint(*this, qoi_index))
1086  );
1087  }
1088  }
1089 
1090 #endif // LIBMESH_ENABLE_DIRICHLET
1091 
1092  STOP_LOG("create_dof_constraints()", "DofMap");
1093 }
void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited

Definition at line 106 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

107 {
108  _enable_print_counter = false;
109  return;
110 }
void libMesh::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 841 of file dof_map.C.

References _end_df, _end_old_df, _first_df, _first_old_df, _n_dfs, _n_old_dfs, _send_list, add_neighbors_to_send_list(), libMesh::Parallel::Communicator::allgather(), libMesh::ParallelObject::comm(), distribute_local_dofs_node_major(), distribute_local_dofs_var_major(), libMesh::DofObject::dof_number(), elem_ptr(), libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), end_dof(), first_dof(), invalidate_dofs(), libMesh::MeshBase::is_prepared(), libMesh::libmesh_assert(), libMesh::MeshTools::libmesh_assert_valid_dof_ids(), mesh, libMesh::DofObject::n_comp(), libMesh::ParallelObject::n_processors(), libMesh::DofObject::n_vars(), node_ptr(), libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), libMesh::on_command_line(), libMesh::ParallelObject::processor_id(), libMesh::DofObject::processor_id(), reinit(), set_nonlocal_dof_objects(), libMesh::START_LOG(), libMesh::STOP_LOG(), and sys_number().

Referenced by libMesh::EquationSystems::allgather(), and libMesh::EquationSystems::reinit().

842 {
843  // This function must be run on all processors at once
844  parallel_object_only();
845 
846  // Log how long it takes to distribute the degrees of freedom
847  START_LOG("distribute_dofs()", "DofMap");
848 
849  libmesh_assert (mesh.is_prepared());
850 
851  const processor_id_type proc_id = this->processor_id();
852  const processor_id_type n_proc = this->n_processors();
853 
854 // libmesh_assert_greater (this->n_variables(), 0);
855  libmesh_assert_less (proc_id, n_proc);
856 
857  // re-init in case the mesh has changed
858  this->reinit(mesh);
859 
860  // By default distribute variables in a
861  // var-major fashion, but allow run-time
862  // specification
863  bool node_major_dofs = libMesh::on_command_line ("--node_major_dofs");
864 
865  // The DOF counter, will be incremented as we encounter
866  // new degrees of freedom
867  dof_id_type next_free_dof = 0;
868 
869  // Clear the send list before we rebuild it
870  _send_list.clear();
871 
872  // Set temporary DOF indices on this processor
873  if (node_major_dofs)
874  this->distribute_local_dofs_node_major (next_free_dof, mesh);
875  else
876  this->distribute_local_dofs_var_major (next_free_dof, mesh);
877 
878  // Get DOF counts on all processors
879  std::vector<dof_id_type> dofs_on_proc(n_proc, 0);
880  this->comm().allgather(next_free_dof, dofs_on_proc);
881 
882  // Resize and fill the _first_df and _end_df arrays
883 #ifdef LIBMESH_ENABLE_AMR
886 #endif
887 
888  _first_df.resize(n_proc);
889  _end_df.resize (n_proc);
890 
891  // Get DOF offsets
892  _first_df[0] = 0;
893  for (processor_id_type i=1; i < n_proc; ++i)
894  _first_df[i] = _end_df[i-1] = _first_df[i-1] + dofs_on_proc[i-1];
895  _end_df[n_proc-1] = _first_df[n_proc-1] + dofs_on_proc[n_proc-1];
896 
897  // Clear all the current DOF indices
898  // (distribute_dofs expects them cleared!)
899  this->invalidate_dofs(mesh);
900 
901  next_free_dof = _first_df[proc_id];
902 
903  // Set permanent DOF indices on this processor
904  if (node_major_dofs)
905  this->distribute_local_dofs_node_major (next_free_dof, mesh);
906  else
907  this->distribute_local_dofs_var_major (next_free_dof, mesh);
908 
909  libmesh_assert_equal_to (next_free_dof, _end_df[proc_id]);
910 
911  //------------------------------------------------------------
912  // At this point, all n_comp and dof_number values on local
913  // DofObjects should be correct, but a ParallelMesh might have
914  // incorrect values on non-local DofObjects. Let's request the
915  // correct values from each other processor.
916 
917  if (this->n_processors() > 1)
918  {
919  this->set_nonlocal_dof_objects(mesh.nodes_begin(),
920  mesh.nodes_end(),
922 
923  this->set_nonlocal_dof_objects(mesh.elements_begin(),
924  mesh.elements_end(),
926  }
927 
928 #ifdef DEBUG
929  {
930  const unsigned int
931  sys_num = this->sys_number();
932 
933  // Processors should all agree on DoF ids
935 
936  // DoF processor ids should match DofObject processor ids
937  MeshBase::const_node_iterator node_it = mesh.nodes_begin();
938  const MeshBase::const_node_iterator node_end = mesh.nodes_end();
939  for ( ; node_it != node_end; ++node_it)
940  {
941  DofObject const * const dofobj = *node_it;
942  const processor_id_type proc_id = dofobj->processor_id();
943 
944  for (unsigned int v=0; v != dofobj->n_vars(sys_num); ++v)
945  for (unsigned int c=0; c != dofobj->n_comp(sys_num,v); ++c)
946  {
947  const dof_id_type dofid = dofobj->dof_number(sys_num,v,c);
948  libmesh_assert_greater_equal (dofid, this->first_dof(proc_id));
949  libmesh_assert_less (dofid, this->end_dof(proc_id));
950  }
951  }
952 
953  MeshBase::const_element_iterator elem_it = mesh.elements_begin();
954  const MeshBase::const_element_iterator elem_end = mesh.elements_end();
955  for ( ; elem_it != elem_end; ++elem_it)
956  {
957  DofObject const * const dofobj = *elem_it;
958  const processor_id_type proc_id = dofobj->processor_id();
959 
960  for (unsigned int v=0; v != dofobj->n_vars(sys_num); ++v)
961  for (unsigned int c=0; c != dofobj->n_comp(sys_num,v); ++c)
962  {
963  const dof_id_type dofid = dofobj->dof_number(sys_num,v,c);
964  libmesh_assert_greater_equal (dofid, this->first_dof(proc_id));
965  libmesh_assert_less (dofid, this->end_dof(proc_id));
966  }
967  }
968  }
969 #endif
970 
971  // Set the total number of degrees of freedom
972 #ifdef LIBMESH_ENABLE_AMR
973  _n_old_dfs = _n_dfs;
974 #endif
975  _n_dfs = _end_df[n_proc-1];
976 
977  STOP_LOG("distribute_dofs()", "DofMap");
978 
979  // Note that in the add_neighbors_to_send_list nodes on processor
980  // boundaries that are shared by multiple elements are added for
981  // each element.
983 
984  // Here we used to clean up that data structure; now System and
985  // EquationSystems call that for us, after we've added constraint
986  // dependencies to the send_list too.
987  // this->sort_send_list ();
988 }
void libMesh::DofMap::distribute_local_dofs_node_major ( dof_id_type 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 991 of file dof_map.C.

References _n_SCALAR_dofs, libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::Variable::active_on_subdomain(), libMesh::FEType::family, libMesh::Elem::get_node(), libMesh::DofObject::invalid_id, libMesh::libmesh_assert(), libMesh::MeshTools::libmesh_assert_valid_procids< Node >(), libMesh::MeshBase::local_nodes_begin(), libMesh::MeshBase::local_nodes_end(), mesh, libMesh::DofObject::n_comp(), libMesh::DofObject::n_comp_group(), n_nodes, libMesh::Elem::n_nodes(), libMesh::ParallelObject::n_processors(), libMesh::DofObject::n_var_groups(), n_variable_groups(), libMesh::VariableGroup::n_variables(), libMesh::FEType::order, libMesh::ParallelObject::processor_id(), libMesh::DofObject::processor_id(), libMeshEnums::SCALAR, libMesh::DofObject::set_vg_dof_base(), libMesh::Elem::subdomain_id(), sys_number(), libMesh::Variable::type(), variable_group(), and libMesh::DofObject::vg_dof_base().

Referenced by distribute_dofs().

993 {
994  const unsigned int sys_num = this->sys_number();
995  const unsigned int n_var_groups = this->n_variable_groups();
996 
997  //-------------------------------------------------------------------------
998  // First count and assign temporary numbers to local dofs
999  MeshBase::element_iterator elem_it = mesh.active_local_elements_begin();
1000  const MeshBase::element_iterator elem_end = mesh.active_local_elements_end();
1001 
1002  for ( ; elem_it != elem_end; ++elem_it)
1003  {
1004  // Only number dofs connected to active
1005  // elements on this processor.
1006  Elem* elem = *elem_it;
1007  const unsigned int n_nodes = elem->n_nodes();
1008 
1009  // First number the nodal DOFS
1010  for (unsigned int n=0; n<n_nodes; n++)
1011  {
1012  Node* node = elem->get_node(n);
1013 
1014  for (unsigned vg=0; vg<n_var_groups; vg++)
1015  {
1016  const VariableGroup &vg_description(this->variable_group(vg));
1017 
1018  if( (vg_description.type().family != SCALAR) &&
1019  (vg_description.active_on_subdomain(elem->subdomain_id())) )
1020  {
1021  // assign dof numbers (all at once) if this is
1022  // our node and if they aren't already there
1023  if ((node->n_comp_group(sys_num,vg) > 0) &&
1024  (node->processor_id() == this->processor_id()) &&
1025  (node->vg_dof_base(sys_num,vg) ==
1027  {
1028  node->set_vg_dof_base(sys_num,
1029  vg,
1030  next_free_dof);
1031  next_free_dof += (vg_description.n_variables()*
1032  node->n_comp_group(sys_num,vg));
1033  //node->debug_buffer();
1034  }
1035  }
1036  }
1037  }
1038 
1039  // Now number the element DOFS
1040  for (unsigned vg=0; vg<n_var_groups; vg++)
1041  {
1042  const VariableGroup &vg_description(this->variable_group(vg));
1043 
1044  if ( (vg_description.type().family != SCALAR) &&
1045  (vg_description.active_on_subdomain(elem->subdomain_id())) )
1046  if (elem->n_comp_group(sys_num,vg) > 0)
1047  {
1048  libmesh_assert_equal_to (elem->vg_dof_base(sys_num,vg),
1050 
1051  elem->set_vg_dof_base(sys_num,
1052  vg,
1053  next_free_dof);
1054 
1055  next_free_dof += (vg_description.n_variables()*
1056  elem->n_comp(sys_num,vg));
1057  }
1058  }
1059  } // done looping over elements
1060 
1061 
1062  // we may have missed assigning DOFs to nodes that we own
1063  // but to which we have no connected elements matching our
1064  // variable restriction criterion. this will happen, for example,
1065  // if variable V is restricted to subdomain S. We may not own
1066  // any elements which live in S, but we may own nodes which are
1067  // *connected* to elements which do. in this scenario these nodes
1068  // will presently have unnumbered DOFs. we need to take care of
1069  // them here since we own them and no other processor will touch them.
1070  {
1071  MeshBase::node_iterator node_it = mesh.local_nodes_begin();
1072  const MeshBase::node_iterator node_end = mesh.local_nodes_end();
1073 
1074  for (; node_it != node_end; ++node_it)
1075  {
1076  Node *node = *node_it;
1077  libmesh_assert(node);
1078 
1079  for (unsigned vg=0; vg<n_var_groups; vg++)
1080  {
1081  const VariableGroup &vg_description(this->variable_group(vg));
1082 
1083  if (node->n_comp_group(sys_num,vg))
1084  if (node->vg_dof_base(sys_num,vg) == DofObject::invalid_id)
1085  {
1086  node->set_vg_dof_base (sys_num,
1087  vg,
1088  next_free_dof);
1089 
1090  next_free_dof += (vg_description.n_variables()*
1091  node->n_comp(sys_num,vg));
1092  }
1093  }
1094  }
1095  }
1096 
1097  // Finally, count up the SCALAR dofs
1098  this->_n_SCALAR_dofs = 0;
1099  for (unsigned vg=0; vg<n_var_groups; vg++)
1100  {
1101  const VariableGroup &vg_description(this->variable_group(vg));
1102 
1103  if( vg_description.type().family == SCALAR )
1104  {
1105  this->_n_SCALAR_dofs += (vg_description.n_variables()*
1106  vg_description.type().order);
1107  continue;
1108  }
1109  }
1110 
1111  // Only increment next_free_dof if we're on the processor
1112  // that holds this SCALAR variable
1113  if ( this->processor_id() == (this->n_processors()-1) )
1114  next_free_dof += _n_SCALAR_dofs;
1115 
1116 #ifdef DEBUG
1117  {
1118  // libMesh::out << "next_free_dof=" << next_free_dof << std::endl
1119  // << "_n_SCALAR_dofs=" << _n_SCALAR_dofs << std::endl;
1120 
1121  // Make sure we didn't miss any nodes
1123 
1124  MeshBase::node_iterator node_it = mesh.local_nodes_begin();
1125  const MeshBase::node_iterator node_end = mesh.local_nodes_end();
1126  for (; node_it != node_end; ++node_it)
1127  {
1128  Node *obj = *node_it;
1129  libmesh_assert(obj);
1130  unsigned int n_var_g = obj->n_var_groups(this->sys_number());
1131  for (unsigned int vg=0; vg != n_var_g; ++vg)
1132  {
1133  unsigned int n_comp_g =
1134  obj->n_comp_group(this->sys_number(), vg);
1135  dof_id_type my_first_dof = n_comp_g ?
1136  obj->vg_dof_base(this->sys_number(), vg) : 0;
1137  libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
1138  }
1139  }
1140  }
1141 #endif // DEBUG
1142 }
void libMesh::DofMap::distribute_local_dofs_var_major ( dof_id_type 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 1146 of file dof_map.C.

References _n_SCALAR_dofs, libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::Variable::active_on_subdomain(), libMesh::FEType::family, libMesh::Elem::get_node(), libMesh::DofObject::invalid_id, libMesh::libmesh_assert(), libMesh::MeshTools::libmesh_assert_valid_procids< Node >(), libMesh::MeshBase::local_nodes_begin(), libMesh::MeshBase::local_nodes_end(), mesh, libMesh::DofObject::n_comp_group(), n_nodes, libMesh::Elem::n_nodes(), libMesh::ParallelObject::n_processors(), libMesh::DofObject::n_var_groups(), n_variable_groups(), libMesh::VariableGroup::n_variables(), libMesh::FEType::order, libMesh::ParallelObject::processor_id(), libMesh::DofObject::processor_id(), libMeshEnums::SCALAR, libMesh::DofObject::set_vg_dof_base(), libMesh::Elem::subdomain_id(), sys_number(), libMesh::Variable::type(), variable_group(), and libMesh::DofObject::vg_dof_base().

Referenced by distribute_dofs().

1148 {
1149  const unsigned int sys_num = this->sys_number();
1150  const unsigned int n_var_groups = this->n_variable_groups();
1151 
1152  //-------------------------------------------------------------------------
1153  // First count and assign temporary numbers to local dofs
1154  for (unsigned vg=0; vg<n_var_groups; vg++)
1155  {
1156  const VariableGroup &vg_description(this->variable_group(vg));
1157 
1158  const unsigned int n_vars_in_group = vg_description.n_variables();
1159 
1160  // Skip the SCALAR dofs
1161  if (vg_description.type().family == SCALAR)
1162  continue;
1163 
1164  MeshBase::element_iterator elem_it = mesh.active_local_elements_begin();
1165  const MeshBase::element_iterator elem_end = mesh.active_local_elements_end();
1166 
1167  for ( ; elem_it != elem_end; ++elem_it)
1168  {
1169  // Only number dofs connected to active
1170  // elements on this processor.
1171  Elem* elem = *elem_it;
1172 
1173  // ... and only variables which are active on
1174  // on this element's subdomain
1175  if (!vg_description.active_on_subdomain(elem->subdomain_id()))
1176  continue;
1177 
1178  const unsigned int n_nodes = elem->n_nodes();
1179 
1180  // First number the nodal DOFS
1181  for (unsigned int n=0; n<n_nodes; n++)
1182  {
1183  Node* node = elem->get_node(n);
1184 
1185  // assign dof numbers (all at once) if this is
1186  // our node and if they aren't already there
1187  if ((node->n_comp_group(sys_num,vg) > 0) &&
1188  (node->processor_id() == this->processor_id()) &&
1189  (node->vg_dof_base(sys_num,vg) ==
1191  {
1192  node->set_vg_dof_base(sys_num,
1193  vg,
1194  next_free_dof);
1195 
1196  next_free_dof += (n_vars_in_group*
1197  node->n_comp_group(sys_num,vg));
1198  }
1199  }
1200 
1201  // Now number the element DOFS
1202  if (elem->n_comp_group(sys_num,vg) > 0)
1203  {
1204  libmesh_assert_equal_to (elem->vg_dof_base(sys_num,vg),
1206 
1207  elem->set_vg_dof_base(sys_num,
1208  vg,
1209  next_free_dof);
1210 
1211  next_free_dof += (n_vars_in_group*
1212  elem->n_comp_group(sys_num,vg));
1213  }
1214  } // end loop on elements
1215 
1216  // we may have missed assigning DOFs to nodes that we own
1217  // but to which we have no connected elements matching our
1218  // variable restriction criterion. this will happen, for example,
1219  // if variable V is restricted to subdomain S. We may not own
1220  // any elements which live in S, but we may own nodes which are
1221  // *connected* to elements which do. in this scenario these nodes
1222  // will presently have unnumbered DOFs. we need to take care of
1223  // them here since we own them and no other processor will touch them.
1224  {
1225  MeshBase::node_iterator node_it = mesh.local_nodes_begin();
1226  const MeshBase::node_iterator node_end = mesh.local_nodes_end();
1227 
1228  for (; node_it != node_end; ++node_it)
1229  {
1230  Node *node = *node_it;
1231  libmesh_assert(node);
1232 
1233  if (node->n_comp_group(sys_num,vg))
1234  if (node->vg_dof_base(sys_num,vg) == DofObject::invalid_id)
1235  {
1236  node->set_vg_dof_base (sys_num,
1237  vg,
1238  next_free_dof);
1239 
1240  next_free_dof += (n_vars_in_group*
1241  node->n_comp_group(sys_num,vg));
1242  }
1243  }
1244  }
1245  } // end loop on variable groups
1246 
1247  // Finally, count up the SCALAR dofs
1248  this->_n_SCALAR_dofs = 0;
1249  for (unsigned vg=0; vg<n_var_groups; vg++)
1250  {
1251  const VariableGroup &vg_description(this->variable_group(vg));
1252 
1253  if( vg_description.type().family == SCALAR )
1254  {
1255  this->_n_SCALAR_dofs += (vg_description.n_variables()*
1256  vg_description.type().order);
1257  continue;
1258  }
1259  }
1260 
1261  // Only increment next_free_dof if we're on the processor
1262  // that holds this SCALAR variable
1263  if ( this->processor_id() == (this->n_processors()-1) )
1264  next_free_dof += _n_SCALAR_dofs;
1265 
1266 #ifdef DEBUG
1267  {
1268  // Make sure we didn't miss any nodes
1270 
1271  MeshBase::node_iterator node_it = mesh.local_nodes_begin();
1272  const MeshBase::node_iterator node_end = mesh.local_nodes_end();
1273  for (; node_it != node_end; ++node_it)
1274  {
1275  Node *obj = *node_it;
1276  libmesh_assert(obj);
1277  unsigned int n_var_g = obj->n_var_groups(this->sys_number());
1278  for (unsigned int vg=0; vg != n_var_g; ++vg)
1279  {
1280  unsigned int n_comp_g =
1281  obj->n_comp_group(this->sys_number(), vg);
1282  dof_id_type my_first_dof = n_comp_g ?
1283  obj->vg_dof_base(this->sys_number(), vg) : 0;
1284  libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
1285  }
1286  }
1287  }
1288 #endif // DEBUG
1289 }
void libMesh::DofMap::dof_indices ( const Elem *const  elem,
std::vector< dof_id_type > &  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 1638 of file dof_map.C.

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

Referenced by libMesh::ExactSolution::_compute_error(), libMesh::UniformRefinementEstimator::_estimate_error(), add_neighbors_to_send_list(), libMesh::HPCoarsenTest::add_projection(), libMesh::EquationSystems::build_discontinuous_solution_vector(), libMesh::EquationSystems::build_solution_vector(), libMesh::System::calculate_norm(), libMesh::FEGenericBase< T >::coarsened_dof_values(), libMesh::FEGenericBase< T >::compute_periodic_constraints(), libMesh::FEGenericBase< T >::compute_proj_constraints(), DMCreateDomainDecomposition_libMesh(), DMCreateFieldDecomposition_libMesh(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::for(), libMesh::EquationSystems::get_solution(), libMesh::MeshFunction::gradient(), libMesh::MeshFunction::hessian(), libMesh::SystemSubsetBySubdomain::init(), libMesh::System::local_dof_indices(), libMesh::DGFEMContext::neighbor_side_fe_reinit(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::SparsityPattern::Build::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::MeshFunction::operator()(), libMesh::ProjectSolution::operator()(), libMesh::BoundaryProjectSolution::operator()(), libMesh::ErrorVector::plot_error(), libMesh::FEMContext::pre_fe_reinit(), libMesh::HPCoarsenTest::select_refinement(), libMesh::EnsightIO::write_scalar_ascii(), and libMesh::EnsightIO::write_vector_ascii().

1641 {
1642  START_LOG("dof_indices()", "DofMap");
1643 
1644  libmesh_assert(elem);
1645 
1646  const unsigned int n_nodes = elem->n_nodes();
1647  const ElemType type = elem->type();
1648  const unsigned int sys_num = this->sys_number();
1649  const unsigned int n_vars = this->n_variables();
1650  const unsigned int dim = elem->dim();
1651 
1652  // Clear the DOF indices vector
1653  di.clear();
1654 
1655 #ifdef DEBUG
1656  // Check that sizes match in DEBUG mode
1657  unsigned int tot_size = 0;
1658 #endif
1659 
1660  // Create a vector to indicate which
1661  // SCALAR variables have been requested
1662  std::vector<unsigned int> SCALAR_var_numbers;
1663  SCALAR_var_numbers.clear();
1664 
1665  // Get the dof numbers
1666  for (unsigned int v=0; v<n_vars; v++)
1667  if ((v == vn) || (vn == libMesh::invalid_uint))
1668  {
1669  if(this->variable(v).type().family == SCALAR)
1670  {
1671  // We asked for this variable, so add it to the vector.
1672  SCALAR_var_numbers.push_back(v);
1673 
1674 #ifdef DEBUG
1675  FEType fe_type = this->variable_type(v);
1676  tot_size += FEInterface::n_dofs(dim,
1677  fe_type,
1678  type);
1679 #endif
1680  }
1681  else
1682  if (this->variable(v).active_on_subdomain(elem->subdomain_id()))
1683  { // Do this for all the variables if one was not specified
1684  // or just for the specified variable
1685 
1686  // Increase the polynomial order on p refined elements
1687  FEType fe_type = this->variable_type(v);
1688  fe_type.order = static_cast<Order>(fe_type.order +
1689  elem->p_level());
1690 
1691  const bool extra_hanging_dofs =
1693 
1694 #ifdef DEBUG
1695  tot_size += FEInterface::n_dofs(dim,
1696  fe_type,
1697  type);
1698 #endif
1699 
1700  // Get the node-based DOF numbers
1701  for (unsigned int n=0; n<n_nodes; n++)
1702  {
1703  const Node* node = elem->get_node(n);
1704 
1705  // There is a potential problem with h refinement. Imagine a
1706  // quad9 that has a linear FE on it. Then, on the hanging side,
1707  // it can falsely identify a DOF at the mid-edge node. This is why
1708  // we call FEInterface instead of node->n_comp() directly.
1709  const unsigned int nc = FEInterface::n_dofs_at_node (dim,
1710  fe_type,
1711  type,
1712  n);
1713 
1714  // If this is a non-vertex on a hanging node with extra
1715  // degrees of freedom, we use the non-vertex dofs (which
1716  // come in reverse order starting from the end, to
1717  // simplify p refinement)
1718  if (extra_hanging_dofs && !elem->is_vertex(n))
1719  {
1720  const int dof_offset = node->n_comp(sys_num,v) - nc;
1721 
1722  // We should never have fewer dofs than necessary on a
1723  // node unless we're getting indices on a parent element,
1724  // and we should never need the indices on such a node
1725  if (dof_offset < 0)
1726  {
1727  libmesh_assert(!elem->active());
1728  di.resize(di.size() + nc, DofObject::invalid_id);
1729  }
1730  else
1731  for (int i=node->n_comp(sys_num,v)-1; i>=dof_offset; i--)
1732  {
1733  libmesh_assert_not_equal_to (node->dof_number(sys_num,v,i),
1735  di.push_back(node->dof_number(sys_num,v,i));
1736  }
1737  }
1738  // If this is a vertex or an element without extra hanging
1739  // dofs, our dofs come in forward order coming from the
1740  // beginning
1741  else
1742  for (unsigned int i=0; i<nc; i++)
1743  {
1744  libmesh_assert_not_equal_to (node->dof_number(sys_num,v,i),
1746  di.push_back(node->dof_number(sys_num,v,i));
1747  }
1748  }
1749 
1750  // If there are any element-based DOF numbers, get them
1751  const unsigned int nc = FEInterface::n_dofs_per_elem(dim,
1752  fe_type,
1753  type);
1754  // We should never have fewer dofs than necessary on an
1755  // element unless we're getting indices on a parent element,
1756  // and we should never need those indices
1757  if (nc != 0)
1758  {
1759  if (elem->n_systems() > sys_num &&
1760  nc <= elem->n_comp(sys_num,v))
1761  {
1762  for (unsigned int i=0; i<nc; i++)
1763  {
1764  libmesh_assert_not_equal_to (elem->dof_number(sys_num,v,i),
1766 
1767  di.push_back(elem->dof_number(sys_num,v,i));
1768  }
1769  }
1770  else
1771  {
1772  libmesh_assert(!elem->active() || fe_type.family == LAGRANGE);
1773  di.resize(di.size() + nc, DofObject::invalid_id);
1774  }
1775  }
1776  }
1777  } // end loop over variables
1778 
1779  // Finally append any SCALAR dofs that we asked for.
1780  std::vector<dof_id_type> di_new;
1781  std::vector<unsigned int>::iterator it = SCALAR_var_numbers.begin();
1782  std::vector<unsigned int>::const_iterator it_end = SCALAR_var_numbers.end();
1783  for( ; it != it_end; ++it)
1784  {
1785  this->SCALAR_dof_indices(di_new,*it);
1786  di.insert( di.end(), di_new.begin(), di_new.end());
1787  }
1788 
1789 #ifdef DEBUG
1790  libmesh_assert_equal_to (tot_size, di.size());
1791 #endif
1792 
1793  STOP_LOG("dof_indices()", "DofMap");
1794 }
DofObject * libMesh::DofMap::elem_ptr ( MeshBase mesh,
dof_id_type  i 
) const
private

An adapter function that returns Elem pointers by index

Definition at line 276 of file dof_map.C.

References libMesh::MeshBase::elem().

Referenced by distribute_dofs().

277 {
278  return mesh.elem(i);
279 }
void libMesh::ReferenceCounter::enable_print_counter_info ( )
staticinherited

Methods to enable/disable the reference counter output from print_info()

Definition at line 100 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

101 {
102  _enable_print_counter = true;
103  return;
104 }
dof_id_type libMesh::DofMap::end_dof ( const processor_id_type  proc) 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 488 of file dof_map.h.

References _end_df.

Referenced by DMCreateDomainDecomposition_libMesh(), DMCreateFieldDecomposition_libMesh(), libMesh::SystemSubsetBySubdomain::init(), libMesh::CondensedEigenSystem::initialize_condensed_dofs(), libMesh::SparsityPattern::Build::join(), libMesh::System::local_dof_indices(), libMesh::SparsityPattern::Build::operator()(), and libMesh::SparsityPattern::Build::parallel_sync().

489  { libmesh_assert_less (proc, _end_df.size()); return _end_df[proc]; }
dof_id_type libMesh::DofMap::end_dof ( ) const
inline

Definition at line 491 of file dof_map.h.

References libMesh::ParallelObject::processor_id().

Referenced by add_neighbors_to_send_list(), all_semilocal_indices(), and distribute_dofs().

492  { return this->end_dof(this->processor_id()); }
dof_id_type libMesh::DofMap::end_old_dof ( const processor_id_type  proc) const
inline

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

Definition at line 499 of file dof_map.h.

References _end_old_df.

500  { libmesh_assert_less (proc, _end_old_df.size()); return _end_old_df[proc]; }
dof_id_type libMesh::DofMap::end_old_dof ( ) const
inline

Definition at line 502 of file dof_map.h.

References libMesh::ParallelObject::processor_id().

503  { return this->end_old_dof(this->processor_id()); }
void libMesh::DofMap::enforce_adjoint_constraints_exactly ( NumericVector< Number > &  v,
unsigned int  q 
) const
inline

Heterogenously constrains the numeric vector v, which represents an adjoint solution defined on the mesh for quantity fo interest q. For homogeneous constraints, use enforce_constraints_exactly instead

Definition at line 1941 of file dof_map_constraints.C.

References libMesh::NumericVector< T >::close(), libMesh::NumericVector< T >::closed(), libMesh::comm, libMesh::err, libMesh::AutoPtr< Tp >::get(), libMeshEnums::GHOSTED, libMesh::libmesh_assert(), libMesh::NumericVector< T >::localize(), libMeshEnums::PARALLEL, libMeshEnums::SERIAL, libMesh::NumericVector< T >::set(), libMesh::NumericVector< T >::size(), libMesh::START_LOG(), libMesh::STOP_LOG(), and libMesh::NumericVector< T >::type().

Referenced by libMesh::ImplicitSystem::adjoint_solve().

1943 {
1944  parallel_object_only();
1945 
1946  if (!this->n_constrained_dofs())
1947  return;
1948 
1949  START_LOG("enforce_adjoint_constraints_exactly()","DofMap");
1950 
1951  NumericVector<Number> *v_local = NULL; // will be initialized below
1952  NumericVector<Number> *v_global = NULL; // will be initialized below
1954  if (v.type() == SERIAL)
1955  {
1956  v_built = NumericVector<Number>::build(this->comm());
1957  v_built->init(this->n_dofs(), this->n_local_dofs(), true, PARALLEL);
1958  v_built->close();
1959 
1960  for (dof_id_type i=v_built->first_local_index();
1961  i<v_built->last_local_index(); i++)
1962  v_built->set(i, v(i));
1963  v_built->close();
1964  v_global = v_built.get();
1965 
1966  v_local = &v;
1967  libmesh_assert (v_local->closed());
1968  }
1969  else if (v.type() == PARALLEL)
1970  {
1971  v_built = NumericVector<Number>::build(this->comm());
1972  v_built->init (v.size(), v.size(), true, SERIAL);
1973  v.localize(*v_built);
1974  v_built->close();
1975  v_local = v_built.get();
1976 
1977  v_global = &v;
1978  }
1979  else if (v.type() == GHOSTED)
1980  {
1981  v_local = &v;
1982  v_global = &v;
1983  }
1984  else // unknown v.type()
1985  {
1986  libMesh::err << "ERROR: Unknown v.type() == " << v.type()
1987  << std::endl;
1988  libmesh_error();
1989  }
1990 
1991  // We should never hit these asserts because we should error-out in
1992  // else clause above. Just to be sure we don't try to use v_local
1993  // and v_global uninitialized...
1994  libmesh_assert(v_local);
1995  libmesh_assert(v_global);
1996 
1997  // Do we have any non_homogeneous constraints?
1998  const AdjointDofConstraintValues::const_iterator
1999  adjoint_constraint_map_it = _adjoint_constraint_values.find(q);
2000  const DofConstraintValueMap *constraint_map =
2001  (adjoint_constraint_map_it == _adjoint_constraint_values.end()) ?
2002  NULL : &adjoint_constraint_map_it->second;
2003 
2004  DofConstraints::const_iterator c_it = _dof_constraints.begin();
2005  const DofConstraints::const_iterator c_end = _dof_constraints.end();
2006 
2007  for ( ; c_it != c_end; ++c_it)
2008  {
2009  dof_id_type constrained_dof = c_it->first;
2010  if (constrained_dof < this->first_dof() ||
2011  constrained_dof >= this->end_dof())
2012  continue;
2013 
2014  const DofConstraintRow constraint_row = c_it->second;
2015 
2016  Number exact_value = 0;
2017  if (constraint_map)
2018  {
2019  const DofConstraintValueMap::const_iterator
2020  adjoint_constraint_it =
2021  constraint_map->find(constrained_dof);
2022  if (adjoint_constraint_it != constraint_map->end())
2023  exact_value = adjoint_constraint_it->second;
2024  }
2025 
2026  for (DofConstraintRow::const_iterator
2027  j=constraint_row.begin(); j != constraint_row.end();
2028  ++j)
2029  exact_value += j->second * (*v_local)(j->first);
2030 
2031  v_global->set(constrained_dof, exact_value);
2032  }
2033 
2034  // If the old vector was serial, we probably need to send our values
2035  // to other processors
2036  if (v.type() == SERIAL)
2037  {
2038 #ifndef NDEBUG
2039  v_global->close();
2040 #endif
2041  v_global->localize (v);
2042  }
2043  v.close();
2044 
2045  STOP_LOG("enforce_adjoint_constraints_exactly()","DofMap");
2046 }
void libMesh::DofMap::enforce_constraints_exactly ( const System system,
NumericVector< Number > *  v = NULL,
bool  homogeneous = false 
) 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

If homogeneous == true, heterogeneous constraints are enforced as if they were homogeneous. This might be appropriate for e.g. a vector representing a difference between two heterogeneously-constrained solutions.

Definition at line 1835 of file dof_map_constraints.C.

References libMesh::NumericVector< T >::close(), libMesh::NumericVector< T >::closed(), libMesh::comm, libMesh::err, libMesh::AutoPtr< Tp >::get(), libMesh::System::get_dof_map(), libMeshEnums::GHOSTED, libMesh::libmesh_assert(), libMesh::NumericVector< T >::localize(), libMeshEnums::PARALLEL, libMeshEnums::SERIAL, libMesh::NumericVector< T >::set(), libMesh::NumericVector< T >::size(), libMesh::System::solution, libMesh::START_LOG(), libMesh::STOP_LOG(), and libMesh::NumericVector< T >::type().

Referenced by libMesh::__libmesh_petsc_diff_solver_jacobian(), libMesh::__libmesh_petsc_diff_solver_residual(), libMesh::__libmesh_petsc_snes_jacobian(), libMesh::__libmesh_petsc_snes_residual(), libMesh::Problem_Interface::computeF(), libMesh::Problem_Interface::computeJacobian(), libMesh::Problem_Interface::computePreconditioner(), DMlibMeshFunction(), DMlibMeshJacobian(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::System::project_vector(), libMesh::ImplicitSystem::sensitivity_solve(), libMesh::NewtonSolver::solve(), libMesh::PetscDiffSolver::solve(), libMesh::ImplicitSystem::weighted_sensitivity_adjoint_solve(), and libMesh::ImplicitSystem::weighted_sensitivity_solve().

1838 {
1839  parallel_object_only();
1840 
1841  if (!this->n_constrained_dofs())
1842  return;
1843 
1844  START_LOG("enforce_constraints_exactly()","DofMap");
1845 
1846  if (!v)
1847  v = system.solution.get();
1848 
1849  NumericVector<Number> *v_local = NULL; // will be initialized below
1850  NumericVector<Number> *v_global = NULL; // will be initialized below
1852  if (v->type() == SERIAL)
1853  {
1854  v_built = NumericVector<Number>::build(this->comm());
1855  v_built->init(this->n_dofs(), this->n_local_dofs(), true, PARALLEL);
1856  v_built->close();
1857 
1858  for (dof_id_type i=v_built->first_local_index();
1859  i<v_built->last_local_index(); i++)
1860  v_built->set(i, (*v)(i));
1861  v_built->close();
1862  v_global = v_built.get();
1863 
1864  v_local = v;
1865  libmesh_assert (v_local->closed());
1866  }
1867  else if (v->type() == PARALLEL)
1868  {
1869  v_built = NumericVector<Number>::build(this->comm());
1870  v_built->init (v->size(), v->size(), true, SERIAL);
1871  v->localize(*v_built);
1872  v_built->close();
1873  v_local = v_built.get();
1874 
1875  v_global = v;
1876  }
1877  else if (v->type() == GHOSTED)
1878  {
1879  v_local = v;
1880  v_global = v;
1881  }
1882  else // unknown v->type()
1883  {
1884  libMesh::err << "ERROR: Unknown v->type() == " << v->type()
1885  << std::endl;
1886  libmesh_error();
1887  }
1888 
1889  // We should never hit these asserts because we should error-out in
1890  // else clause above. Just to be sure we don't try to use v_local
1891  // and v_global uninitialized...
1892  libmesh_assert(v_local);
1893  libmesh_assert(v_global);
1894  libmesh_assert_equal_to (this, &(system.get_dof_map()));
1895 
1896  DofConstraints::const_iterator c_it = _dof_constraints.begin();
1897  const DofConstraints::const_iterator c_end = _dof_constraints.end();
1898 
1899  for ( ; c_it != c_end; ++c_it)
1900  {
1901  dof_id_type constrained_dof = c_it->first;
1902  if (constrained_dof < this->first_dof() ||
1903  constrained_dof >= this->end_dof())
1904  continue;
1905 
1906  const DofConstraintRow constraint_row = c_it->second;
1907 
1908  Number exact_value = 0;
1909  if (!homogeneous)
1910  {
1911  DofConstraintValueMap::const_iterator rhsit =
1912  _primal_constraint_values.find(constrained_dof);
1913  if (rhsit != _primal_constraint_values.end())
1914  exact_value = rhsit->second;
1915  }
1916  for (DofConstraintRow::const_iterator
1917  j=constraint_row.begin(); j != constraint_row.end();
1918  ++j)
1919  exact_value += j->second * (*v_local)(j->first);
1920 
1921  v_global->set(constrained_dof, exact_value);
1922  }
1923 
1924  // If the old vector was serial, we probably need to send our values
1925  // to other processors
1926  if (v->type() == SERIAL)
1927  {
1928 #ifndef NDEBUG
1929  v_global->close();
1930 #endif
1931  v_global->localize (*v);
1932  }
1933  v->close();
1934 
1935  STOP_LOG("enforce_constraints_exactly()","DofMap");
1936 }
void libMesh::DofMap::extract_local_vector ( const NumericVector< Number > &  Ug,
const std::vector< dof_id_type > &  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 $ Ue[i] = Ug[dof_indices[i]] $

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 1552 of file dof_map.C.

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

1555 {
1556 #ifdef LIBMESH_ENABLE_AMR
1557 
1558  // Trivial mapping
1559  libmesh_assert_equal_to (dof_indices_in.size(), Ue.size());
1560  bool has_constrained_dofs = false;
1561 
1562  for (unsigned int il=0;
1563  il != libmesh_cast_int<unsigned int>(dof_indices_in.size()); il++)
1564  {
1565  const dof_id_type ig = dof_indices_in[il];
1566 
1567  if (this->is_constrained_dof (ig)) has_constrained_dofs = true;
1568 
1569  libmesh_assert_less (ig, Ug.size());
1570 
1571  Ue.el(il) = Ug(ig);
1572  }
1573 
1574  // If the element has any constrained DOFs then we need
1575  // to account for them in the mapping. This will handle
1576  // the case that the input vector is not constrained.
1577  if (has_constrained_dofs)
1578  {
1579  // Copy the input DOF indices.
1580  std::vector<dof_id_type> constrained_dof_indices(dof_indices_in);
1581 
1582  DenseMatrix<Number> C;
1583  DenseVector<Number> H;
1584 
1585  this->build_constraint_matrix_and_vector (C, H, constrained_dof_indices);
1586 
1587  libmesh_assert_equal_to (dof_indices_in.size(), C.m());
1588  libmesh_assert_equal_to (constrained_dof_indices.size(), C.n());
1589 
1590  // zero-out Ue
1591  Ue.zero();
1592 
1593  // compute Ue = C Ug, with proper mapping.
1594  const unsigned int n_original_dofs =
1595  libmesh_cast_int<unsigned int>(dof_indices_in.size());
1596  for (unsigned int i=0; i != n_original_dofs; i++)
1597  {
1598  Ue.el(i) = H(i);
1599 
1600  const unsigned int n_constrained =
1601  libmesh_cast_int<unsigned int>(constrained_dof_indices.size());
1602  for (unsigned int j=0; j<n_constrained; j++)
1603  {
1604  const dof_id_type jg = constrained_dof_indices[j];
1605 
1606 // If Ug is a serial or ghosted vector, then this assert is
1607 // overzealous. If Ug is a parallel vector, then this assert
1608 // is redundant.
1609 // libmesh_assert ((jg >= Ug.first_local_index()) &&
1610 // (jg < Ug.last_local_index()));
1611 
1612  Ue.el(i) += C(i,j)*Ug(jg);
1613  }
1614  }
1615  }
1616 
1617 #else
1618 
1619  // Trivial mapping
1620 
1621  const unsigned int n_original_dofs =
1622  libmesh_cast_int<unsigned int>(dof_indices_in.size());
1623 
1624  libmesh_assert_equal_to (n_original_dofs, Ue.size());
1625 
1626  for (unsigned int il=0; il<n_original_dofs; il++)
1627  {
1628  const dof_id_type ig = dof_indices_in[il];
1629 
1630  libmesh_assert ((ig >= Ug.first_local_index()) && (ig < Ug.last_local_index()));
1631 
1632  Ue.el(il) = Ug(ig);
1633  }
1634 
1635 #endif
1636 }
void libMesh::DofMap::find_connected_dof_objects ( std::vector< const DofObject * > &  objs) const
private

Finds all the DofObjects associated with the set in objs. This will account for off-element couplings via hanging nodes.

void libMesh::DofMap::find_connected_dofs ( std::vector< dof_id_type > &  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 2128 of file dof_map.C.

References _dof_constraints, is_constrained_dof(), and libMesh::libmesh_assert().

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

2129 {
2130  typedef std::set<dof_id_type> RCSet;
2131 
2132  // First insert the DOFS we already depend on into the set.
2133  RCSet dof_set (elem_dofs.begin(), elem_dofs.end());
2134 
2135  bool done = true;
2136 
2137  // Next insert any dofs those might be constrained in terms
2138  // of. Note that in this case we may not be done: Those may
2139  // in turn depend on others. So, we need to repeat this process
2140  // in that case until the system depends only on unconstrained
2141  // degrees of freedom.
2142  for (unsigned int i=0; i<elem_dofs.size(); i++)
2143  if (this->is_constrained_dof(elem_dofs[i]))
2144  {
2145  // If the DOF is constrained
2146  DofConstraints::const_iterator
2147  pos = _dof_constraints.find(elem_dofs[i]);
2148 
2149  libmesh_assert (pos != _dof_constraints.end());
2150 
2151  const DofConstraintRow& constraint_row = pos->second;
2152 
2153 // adaptive p refinement currently gives us lots of empty constraint
2154 // rows - we should optimize those DoFs away in the future. [RHS]
2155 // libmesh_assert (!constraint_row.empty());
2156 
2157  DofConstraintRow::const_iterator it = constraint_row.begin();
2158  DofConstraintRow::const_iterator it_end = constraint_row.end();
2159 
2160 
2161  // Add the DOFs this dof is constrained in terms of.
2162  // note that these dofs might also be constrained, so
2163  // we will need to call this function recursively.
2164  for ( ; it != it_end; ++it)
2165  if (!dof_set.count (it->first))
2166  {
2167  dof_set.insert (it->first);
2168  done = false;
2169  }
2170  }
2171 
2172 
2173  // If not done then we need to do more work
2174  // (obviously :-) )!
2175  if (!done)
2176  {
2177  // Fill the vector with the contents of the set
2178  elem_dofs.clear();
2179  elem_dofs.insert (elem_dofs.end(),
2180  dof_set.begin(), dof_set.end());
2181 
2182 
2183  // May need to do this recursively. It is possible
2184  // that we just replaced a constrained DOF with another
2185  // constrained DOF.
2186  this->find_connected_dofs (elem_dofs);
2187 
2188  } // end if (!done)
2189 }
dof_id_type libMesh::DofMap::first_dof ( const processor_id_type  proc) const
inline
dof_id_type libMesh::DofMap::first_dof ( ) const
inline

Definition at line 455 of file dof_map.h.

References libMesh::ParallelObject::processor_id().

Referenced by add_neighbors_to_send_list(), all_semilocal_indices(), and distribute_dofs().

456  { return this->first_dof(this->processor_id()); }
dof_id_type libMesh::DofMap::first_old_dof ( const processor_id_type  proc) const
inline

Returns the first old dof index that is local to partition proc.

Definition at line 462 of file dof_map.h.

References _first_old_df.

463  { libmesh_assert_less (proc, _first_old_df.size()); return _first_old_df[proc]; }
dof_id_type libMesh::DofMap::first_old_dof ( ) const
inline

Definition at line 465 of file dof_map.h.

References libMesh::ParallelObject::processor_id().

466  { return this->first_old_dof(this->processor_id()); }
const DirichletBoundaries * libMesh::DofMap::get_adjoint_dirichlet_boundaries ( unsigned int  q) const

Definition at line 3583 of file dof_map_constraints.C.

References libMesh::libmesh_assert_greater().

3584  {
3587  }
DirichletBoundaries * libMesh::DofMap::get_adjoint_dirichlet_boundaries ( unsigned int  q)

Definition at line 3591 of file dof_map_constraints.C.

3592  {
3593  std::size_t old_size = _adjoint_dirichlet_boundaries.size();
3594  for (std::size_t i = old_size; i <= q; ++i)
3596 
3598  }
const DirichletBoundaries* libMesh::DofMap::get_dirichlet_boundaries ( ) const
inline

Definition at line 942 of file dof_map.h.

References _dirichlet_boundaries.

943  {
944  return _dirichlet_boundaries;
945  }
DirichletBoundaries* libMesh::DofMap::get_dirichlet_boundaries ( )
inline

Definition at line 947 of file dof_map.h.

References _dirichlet_boundaries.

948  {
949  return _dirichlet_boundaries;
950  }
std::string libMesh::ReferenceCounter::get_info ( )
staticinherited

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

References libMesh::ReferenceCounter::_counts, and libMesh::Quality::name().

Referenced by libMesh::ReferenceCounter::print_info().

48 {
49 #if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
50 
51  std::ostringstream oss;
52 
53  oss << '\n'
54  << " ---------------------------------------------------------------------------- \n"
55  << "| Reference count information |\n"
56  << " ---------------------------------------------------------------------------- \n";
57 
58  for (Counts::iterator it = _counts.begin();
59  it != _counts.end(); ++it)
60  {
61  const std::string name(it->first);
62  const unsigned int creations = it->second.first;
63  const unsigned int destructions = it->second.second;
64 
65  oss << "| " << name << " reference count information:\n"
66  << "| Creations: " << creations << '\n'
67  << "| Destructions: " << destructions << '\n';
68  }
69 
70  oss << " ---------------------------------------------------------------------------- \n";
71 
72  return oss.str();
73 
74 #else
75 
76  return "";
77 
78 #endif
79 }
std::string libMesh::DofMap::get_info ( ) const

Gets summary info about the sparsity bandwidth and constraints.

Definition at line 2831 of file dof_map.C.

References libMesh::comm, end, libMesh::libmesh_assert(), std::max(), libMesh::processor_id(), libMesh::DofObject::processor_id(), and libMesh::TypeVector< T >::size().

Referenced by libMesh::System::get_info().

2832 {
2833  std::ostringstream os;
2834 
2835  // If we didn't calculate the exact sparsity pattern, the threaded
2836  // sparsity pattern assembly may have just given us an upper bound
2837  // on sparsity.
2838  const char* may_equal = " <= ";
2839 
2840  // If we calculated the exact sparsity pattern, then we can report
2841  // exact bandwidth figures:
2842  std::vector<SparseMatrix<Number>* >::const_iterator
2843  pos = _matrices.begin(),
2844  end = _matrices.end();
2845 
2846  for (; pos != end; ++pos)
2847  if ((*pos)->need_full_sparsity_pattern())
2848  may_equal = " = ";
2849 
2850  dof_id_type max_n_nz = 0, max_n_oz = 0;
2851  long double avg_n_nz = 0, avg_n_oz = 0;
2852 
2853  if (_n_nz)
2854  {
2855  for (std::size_t i = 0; i != _n_nz->size(); ++i)
2856  {
2857  max_n_nz = std::max(max_n_nz, (*_n_nz)[i]);
2858  avg_n_nz += (*_n_nz)[i];
2859  }
2860 
2861  std::size_t n_nz_size = _n_nz->size();
2862 
2863  this->comm().max(max_n_nz);
2864  this->comm().sum(avg_n_nz);
2865  this->comm().sum(n_nz_size);
2866 
2867  avg_n_nz /= std::max(n_nz_size,std::size_t(1));
2868 
2870 
2871  for (std::size_t i = 0; i != (*_n_oz).size(); ++i)
2872  {
2873  max_n_oz = std::max(max_n_oz, (*_n_oz)[i]);
2874  avg_n_oz += (*_n_oz)[i];
2875  }
2876 
2877  std::size_t n_oz_size = _n_oz->size();
2878 
2879  this->comm().max(max_n_oz);
2880  this->comm().sum(avg_n_oz);
2881  this->comm().sum(n_oz_size);
2882 
2883  avg_n_oz /= std::max(n_oz_size,std::size_t(1));
2884  }
2885 
2886  os << " DofMap Sparsity\n Average On-Processor Bandwidth"
2887  << may_equal << avg_n_nz << '\n';
2888 
2889  os << " Average Off-Processor Bandwidth"
2890  << may_equal << avg_n_oz << '\n';
2891 
2892  os << " Maximum On-Processor Bandwidth"
2893  << may_equal << max_n_nz << '\n';
2894 
2895  os << " Maximum Off-Processor Bandwidth"
2896  << may_equal << max_n_oz << std::endl;
2897 
2898 #ifdef LIBMESH_ENABLE_CONSTRAINTS
2899 
2900  std::size_t n_constraints = 0, max_constraint_length = 0,
2901  n_rhss = 0;
2902  long double avg_constraint_length = 0.;
2903 
2904  for (DofConstraints::const_iterator it=_dof_constraints.begin();
2905  it != _dof_constraints.end(); ++it)
2906  {
2907  // Only count local constraints, then sum later
2908  const dof_id_type constrained_dof = it->first;
2909  if (constrained_dof < this->first_dof() ||
2910  constrained_dof >= this->end_dof())
2911  continue;
2912 
2913  const DofConstraintRow& row = it->second;
2914  std::size_t rowsize = row.size();
2915 
2916  max_constraint_length = std::max(max_constraint_length,
2917  rowsize);
2918  avg_constraint_length += rowsize;
2919  n_constraints++;
2920 
2921  if (_primal_constraint_values.count(constrained_dof))
2922  n_rhss++;
2923  }
2924 
2925  this->comm().sum(n_constraints);
2926  this->comm().sum(n_rhss);
2927  this->comm().sum(avg_constraint_length);
2928  this->comm().max(max_constraint_length);
2929 
2930  os << " DofMap Constraints\n Number of DoF Constraints = "
2931  << n_constraints;
2932  if (n_rhss)
2933  os << '\n'
2934  << " Number of Heterogenous Constraints= " << n_rhss;
2935  if (n_constraints)
2936  {
2937  avg_constraint_length /= n_constraints;
2938 
2939  os << '\n'
2940  << " Average DoF Constraint Length= " << avg_constraint_length;
2941  }
2942 
2943 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2944  std::size_t n_node_constraints = 0, max_node_constraint_length = 0,
2945  n_node_rhss = 0;
2946  long double avg_node_constraint_length = 0.;
2947 
2948  for (NodeConstraints::const_iterator it=_node_constraints.begin();
2949  it != _node_constraints.end(); ++it)
2950  {
2951  // Only count local constraints, then sum later
2952  const Node *node = it->first;
2953  if (node->processor_id() != this->processor_id())
2954  continue;
2955 
2956  const NodeConstraintRow& row = it->second.first;
2957  std::size_t rowsize = row.size();
2958 
2959  max_node_constraint_length = std::max(max_node_constraint_length,
2960  rowsize);
2961  avg_node_constraint_length += rowsize;
2962  n_node_constraints++;
2963 
2964  if (it->second.second != Point(0))
2965  n_node_rhss++;
2966  }
2967 
2968  this->comm().sum(n_node_constraints);
2969  this->comm().sum(n_node_rhss);
2970  this->comm().sum(avg_node_constraint_length);
2971  this->comm().max(max_node_constraint_length);
2972 
2973  os << "\n Number of Node Constraints = " << n_node_constraints;
2974  if (n_node_rhss)
2975  os << '\n'
2976  << " Number of Heterogenous Node Constraints= " << n_node_rhss;
2977  if (n_node_constraints)
2978  {
2979  avg_node_constraint_length /= n_node_constraints;
2980  os << "\n Maximum Node Constraint Length= " << max_node_constraint_length
2981  << '\n'
2982  << " Average Node Constraint Length= " << avg_node_constraint_length;
2983  }
2984 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2985 
2986  os << std::endl;
2987 
2988 #endif // LIBMESH_ENABLE_CONSTRAINTS
2989 
2990  return os.str();
2991 }
std::string libMesh::DofMap::get_local_constraints ( bool  print_nonlocal = false) const

Gets a string reporting all DoF and Node constraints local to this processor. If print_nonlocal is true, then nonlocal constraints which are locally known are included.

Definition at line 1190 of file dof_map_constraints.C.

References libMesh::for(), libMesh::DofObject::id(), libMesh::processor_id(), and libMesh::DofObject::processor_id().

1191 {
1192  std::ostringstream os;
1193 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1194  if (print_nonlocal)
1195  os << "All ";
1196  else
1197  os << "Local ";
1198 
1199  os << "Node Constraints:"
1200  << std::endl;
1201 
1202  for (NodeConstraints::const_iterator it=_node_constraints.begin();
1203  it != _node_constraints.end(); ++it)
1204  {
1205  const Node *node = it->first;
1206 
1207  // Skip non-local nodes if requested
1208  if (!print_nonlocal &&
1209  node->processor_id() != this->processor_id())
1210  continue;
1211 
1212  const NodeConstraintRow& row = it->second.first;
1213  const Point& offset = it->second.second;
1214 
1215  os << "Constraints for Node id " << node->id()
1216  << ": \t";
1217 
1218  for (NodeConstraintRow::const_iterator pos=row.begin();
1219  pos != row.end(); ++pos)
1220  os << " (" << pos->first->id() << ","
1221  << pos->second << ")\t";
1222 
1223  os << "rhs: " << offset;
1224 
1225  os << std::endl;
1226  }
1227 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
1228 
1229  if (print_nonlocal)
1230  os << "All ";
1231  else
1232  os << "Local ";
1233 
1234  os << "DoF Constraints:"
1235  << std::endl;
1236 
1237  for (DofConstraints::const_iterator it=_dof_constraints.begin();
1238  it != _dof_constraints.end(); ++it)
1239  {
1240  const dof_id_type i = it->first;
1241 
1242  // Skip non-local dofs if requested
1243  if (!print_nonlocal &&
1244  ((i < this->first_dof()) ||
1245  (i >= this->end_dof())))
1246  continue;
1247 
1248  const DofConstraintRow& row = it->second;
1249  DofConstraintValueMap::const_iterator rhsit =
1250  _primal_constraint_values.find(i);
1251  const Number rhs = (rhsit == _primal_constraint_values.end()) ?
1252  0 : rhsit->second;
1253 
1254  os << "Constraints for DoF " << i
1255  << ": \t";
1256 
1257  for (DofConstraintRow::const_iterator pos=row.begin();
1258  pos != row.end(); ++pos)
1259  os << " (" << pos->first << ","
1260  << pos->second << ")\t";
1261 
1262  os << "rhs: " << rhs;
1263 
1264  os << std::endl;
1265  }
1266 
1267  for (unsigned int qoi_index = 0;
1268  qoi_index != _adjoint_dirichlet_boundaries.size();
1269  ++qoi_index)
1270  {
1271  os << "Adjoint " << qoi_index << " DoF rhs values:"
1272  << std::endl;
1273 
1274  AdjointDofConstraintValues::const_iterator adjoint_map_it =
1275  _adjoint_constraint_values.find(qoi_index);
1276 
1277  if (adjoint_map_it != _adjoint_constraint_values.end())
1278  for (DofConstraintValueMap::const_iterator
1279  it=adjoint_map_it->second.begin();
1280  it != adjoint_map_it->second.end(); ++it)
1281  {
1282  const dof_id_type i = it->first;
1283 
1284  // Skip non-local dofs if requested
1285  if (!print_nonlocal &&
1286  ((i < this->first_dof()) ||
1287  (i >= this->end_dof())))
1288  continue;
1289 
1290  const Number rhs = it->second;
1291 
1292  os << "RHS for DoF " << i
1293  << ": " << rhs;
1294 
1295  os << std::endl;
1296  }
1297  }
1298 
1299  return os.str();
1300 }
const std::vector<dof_id_type>& libMesh::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 324 of file dof_map.h.

References _n_nz, and libMesh::libmesh_assert().

324  {
326  return *_n_nz;
327  }
const std::vector<dof_id_type>& libMesh::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 335 of file dof_map.h.

References _n_oz, and libMesh::libmesh_assert().

335  {
337  return *_n_oz;
338  }
PeriodicBoundaries* libMesh::DofMap::get_periodic_boundaries ( )
inline

Definition at line 904 of file dof_map.h.

References _periodic_boundaries.

905  {
906  return _periodic_boundaries;
907  }
const std::vector<dof_id_type>& libMesh::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 316 of file dof_map.h.

References _send_list.

Referenced by libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::UnsteadySolver::adjoint_advance_timestep(), libMesh::UnsteadySolver::advance_timestep(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::UnsteadySolver::init_data(), libMesh::System::re_update(), libMesh::UnsteadySolver::reinit(), and libMesh::UnsteadySolver::retrieve_timestep().

316 { return _send_list; }
bool libMesh::DofMap::has_adjoint_dirichlet_boundaries ( unsigned int  q) const

Definition at line 3573 of file dof_map_constraints.C.

Referenced by libMesh::AdjointRefinementEstimator::estimate_error(), and libMesh::ImplicitSystem::weighted_sensitivity_adjoint_solve().

3574  {
3575  if (_adjoint_dirichlet_boundaries.size() > q)
3576  return true;
3577 
3578  return false;
3579  }
bool libMesh::DofMap::has_blocked_representation ( ) const
inline
Returns
true if the variables are capable of being stored in a blocked form. Presently, this means that there can only be one variable group, and that the group has more than one variable.

Definition at line 403 of file dof_map.h.

References n_variable_groups(), and n_variables().

Referenced by block_size().

404  {
405 #ifdef LIBMESH_ENABLE_BLOCKED_STORAGE
406  return ((this->n_variable_groups() == 1) && (this->n_variables() > 1));
407 #else
408  return false;
409 #endif
410  }
bool libMesh::DofMap::has_heterogenous_adjoint_constraint ( const unsigned int  qoi_num,
const dof_id_type  dof 
) const
inline
Returns
true if the degree of freedom dof has a heterogenous constraint for adjoint solution qoi_num, false otherwise.

Definition at line 1441 of file dof_map.h.

References _adjoint_constraint_values.

1443 {
1444  AdjointDofConstraintValues::const_iterator it =
1445  _adjoint_constraint_values.find(qoi_num);
1446  if (it != _adjoint_constraint_values.end() &&
1447  it->second.count(dof))
1448  return true;
1449 
1450  return false;
1451 }
void libMesh::DofMap::heterogenously_constrain_element_matrix_and_vector ( DenseMatrix< Number > &  matrix,
DenseVector< Number > &  rhs,
std::vector< dof_id_type > &  elem_dofs,
bool  asymmetric_constraint_rows = true,
int  qoi_index = -1 
) const

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.

The heterogenous version of this method creates linear systems in which heterogenously constrained degrees of freedom will solve to their correct offset values, as would be appropriate for finding a solution to a linear problem in a single algebraic solve. The non-heterogenous version of this method creates linear systems in which even heterogenously constrained degrees of freedom are solved without offset values taken into account, as would be appropriate for finding linearized updates to a solution in which heterogenous constraints are already satisfied.

By default, the constraints for the primal solution of this system are used. If a non-negative qoi_index is passed in, then the constraints for the corresponding adjoint solution are used instead.

Definition at line 1458 of file dof_map_constraints.C.

References dof_id, libMesh::DenseMatrix< T >::left_multiply_transpose(), libMesh::libmesh_assert(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), libMesh::DenseMatrix< T >::right_multiply(), libMesh::DenseVector< T >::size(), libMesh::START_LOG(), libMesh::STOP_LOG(), libMesh::DenseMatrix< T >::vector_mult(), and libMesh::DenseMatrix< T >::vector_mult_transpose().

1463 {
1464  libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
1465  libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
1466  libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
1467 
1468  // check for easy return
1469  if (this->_dof_constraints.empty())
1470  return;
1471 
1472  // The constrained matrix is built up as C^T K C.
1473  // The constrained RHS is built up as C^T (F - K H)
1476 
1477  this->build_constraint_matrix_and_vector (C, H, elem_dofs, qoi_index);
1478 
1479  START_LOG("hetero_cnstrn_elem_mat_vec()", "DofMap");
1480 
1481  // It is possible that the matrix is not constrained at all.
1482  if ((C.m() == matrix.m()) &&
1483  (C.n() == elem_dofs.size())) // It the matrix is constrained
1484  {
1485  // We may have rhs values to use later
1486  const DofConstraintValueMap *rhs_values = NULL;
1487  if (qoi_index < 0)
1488  rhs_values = &_primal_constraint_values;
1489  else
1490  {
1491  const AdjointDofConstraintValues::const_iterator
1492  it = _adjoint_constraint_values.find(qoi_index);
1493  if (it != _adjoint_constraint_values.end())
1494  rhs_values = &it->second;
1495  }
1496 
1497  // Compute matrix/vector product K H
1499  matrix.vector_mult(KH, H);
1500 
1501  // Compute the matrix-vector product C^T (F - KH)
1502  DenseVector<Number> F_minus_KH(rhs);
1503  F_minus_KH -= KH;
1504  C.vector_mult_transpose(rhs, F_minus_KH);
1505 
1506  // Compute the matrix-matrix-matrix product C^T K C
1507  matrix.left_multiply_transpose (C);
1508  matrix.right_multiply (C);
1509 
1510  libmesh_assert_equal_to (matrix.m(), matrix.n());
1511  libmesh_assert_equal_to (matrix.m(), elem_dofs.size());
1512  libmesh_assert_equal_to (matrix.n(), elem_dofs.size());
1513 
1514  for (unsigned int i=0; i<elem_dofs.size(); i++)
1515  {
1516  const dof_id_type dof_id = elem_dofs[i];
1517 
1518  if (this->is_constrained_dof(dof_id))
1519  {
1520  for (unsigned int j=0; j<matrix.n(); j++)
1521  matrix(i,j) = 0.;
1522 
1523  // If the DOF is constrained
1524  matrix(i,i) = 1.;
1525 
1526  // This will put a nonsymmetric entry in the constraint
1527  // row to ensure that the linear system produces the
1528  // correct value for the constrained DOF.
1529  if (asymmetric_constraint_rows)
1530  {
1531  DofConstraints::const_iterator
1532  pos = _dof_constraints.find(dof_id);
1533 
1534  libmesh_assert (pos != _dof_constraints.end());
1535 
1536  const DofConstraintRow& constraint_row = pos->second;
1537 
1538  for (DofConstraintRow::const_iterator
1539  it=constraint_row.begin(); it != constraint_row.end();
1540  ++it)
1541  for (unsigned int j=0; j<elem_dofs.size(); j++)
1542  if (elem_dofs[j] == it->first)
1543  matrix(i,j) = -it->second;
1544 
1545  if (rhs_values)
1546  {
1547  const DofConstraintValueMap::const_iterator valpos =
1548  rhs_values->find(dof_id);
1549 
1550  rhs(i) = (valpos == rhs_values->end()) ?
1551  0 : valpos->second;
1552  }
1553  }
1554  else
1555  rhs(i) = 0.;
1556  }
1557  }
1558 
1559  } // end if is constrained...
1560 
1561  STOP_LOG("hetero_cnstrn_elem_mat_vec()", "DofMap");
1562 }
void libMesh::DofMap::heterogenously_constrain_element_vector ( const DenseMatrix< Number > &  matrix,
DenseVector< Number > &  rhs,
std::vector< dof_id_type > &  elem_dofs,
bool  asymmetric_constraint_rows = true,
int  qoi_index = -1 
) const

Constrains the element vector. This method requires the element matrix to be square and not-yet-constrained, in which case the elem_dofs correspond to the global DOF indices of both the rows and columns of the element matrix.

The heterogenous version of this method creates linear systems in which heterogenously constrained degrees of freedom will solve to their correct offset values, as would be appropriate for finding a solution to a linear problem in a single algebraic solve. The non-heterogenous version of this method creates linear systems in which even heterogenously constrained degrees of freedom are solved without offset values taken into account, as would be appropriate for finding linearized updates to a solution in which heterogenous constraints are already satisfied.

By default, the constraints for the primal solution of this system are used. If a non-negative qoi_index is passed in, then the constraints for the corresponding adjoint solution are used instead.

Definition at line 1567 of file dof_map_constraints.C.

References dof_id, libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), libMesh::DenseVector< T >::size(), libMesh::START_LOG(), libMesh::STOP_LOG(), libMesh::DenseMatrix< T >::vector_mult(), and libMesh::DenseMatrix< T >::vector_mult_transpose().

1572 {
1573  libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
1574  libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
1575  libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
1576 
1577  // check for easy return
1578  if (this->_dof_constraints.empty())
1579  return;
1580 
1581  // The constrained matrix is built up as C^T K C.
1582  // The constrained RHS is built up as C^T (F - K H)
1585 
1586  this->build_constraint_matrix_and_vector (C, H, elem_dofs, qoi_index);
1587 
1588  START_LOG("hetero_cnstrn_elem_vec()", "DofMap");
1589 
1590  // It is possible that the matrix is not constrained at all.
1591  if ((C.m() == matrix.m()) &&
1592  (C.n() == elem_dofs.size())) // It the matrix is constrained
1593  {
1594  // We may have rhs values to use later
1595  const DofConstraintValueMap *rhs_values = NULL;
1596  if (qoi_index < 0)
1597  rhs_values = &_primal_constraint_values;
1598  else
1599  {
1600  const AdjointDofConstraintValues::const_iterator
1601  it = _adjoint_constraint_values.find(qoi_index);
1602  if (it != _adjoint_constraint_values.end())
1603  rhs_values = &it->second;
1604  }
1605 
1606  // Compute matrix/vector product K H
1608  matrix.vector_mult(KH, H);
1609 
1610  // Compute the matrix-vector product C^T (F - KH)
1611  DenseVector<Number> F_minus_KH(rhs);
1612  F_minus_KH -= KH;
1613  C.vector_mult_transpose(rhs, F_minus_KH);
1614 
1615  for (unsigned int i=0; i<elem_dofs.size(); i++)
1616  {
1617  const dof_id_type dof_id = elem_dofs[i];
1618 
1619  if (this->is_constrained_dof(dof_id))
1620  {
1621  // This will put a nonsymmetric entry in the constraint
1622  // row to ensure that the linear system produces the
1623  // correct value for the constrained DOF.
1624  if (asymmetric_constraint_rows && rhs_values)
1625  {
1626  const DofConstraintValueMap::const_iterator valpos =
1627  rhs_values->find(dof_id);
1628 
1629  rhs(i) = (valpos == rhs_values->end()) ?
1630  0 : valpos->second;
1631  }
1632  else
1633  rhs(i) = 0.;
1634  }
1635  }
1636 
1637  } // end if is constrained...
1638 
1639  STOP_LOG("hetero_cnstrn_elem_vec()", "DofMap");
1640 }
void libMesh::ReferenceCounter::increment_constructor_count ( const std::string &  name)
inlineprotectedinherited

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

Definition at line 163 of file reference_counter.h.

References libMesh::ReferenceCounter::_counts, libMesh::Quality::name(), and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::ReferenceCountedObject().

164 {
165  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
166  std::pair<unsigned int, unsigned int>& p = _counts[name];
167 
168  p.first++;
169 }
void libMesh::ReferenceCounter::increment_destructor_count ( const std::string &  name)
inlineprotectedinherited

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

Definition at line 176 of file reference_counter.h.

References libMesh::ReferenceCounter::_counts, libMesh::Quality::name(), and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::~ReferenceCountedObject().

177 {
178  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
179  std::pair<unsigned int, unsigned int>& p = _counts[name];
180 
181  p.second++;
182 }
void libMesh::DofMap::invalidate_dofs ( MeshBase mesh) const
private

Invalidates all active DofObject dofs for this system

Definition at line 787 of file dof_map.C.

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

Referenced by distribute_dofs(), and reinit().

788 {
789  const unsigned int sys_num = this->sys_number();
790 
791  // All the nodes
792  MeshBase::node_iterator node_it = mesh.nodes_begin();
793  const MeshBase::node_iterator node_end = mesh.nodes_end();
794 
795  for ( ; node_it != node_end; ++node_it)
796  (*node_it)->invalidate_dofs(sys_num);
797 
798  // All the elements
799  MeshBase::element_iterator elem_it = mesh.active_elements_begin();
800  const MeshBase::element_iterator elem_end = mesh.active_elements_end();
801 
802  for ( ; elem_it != elem_end; ++elem_it)
803  (*elem_it)->invalidate_dofs(sys_num);
804 }
bool libMesh::DofMap::is_attached ( SparseMatrix< Number > &  matrix)

Matrices should not be attached more than once. We can test for an already-attached matrix if necessary using is_attached

Definition at line 261 of file dof_map.C.

References _matrices.

Referenced by libMesh::ImplicitSystem::init_matrices().

262 {
263  return (std::find(_matrices.begin(), _matrices.end(),
264  &matrix) != _matrices.end());
265 }
bool libMesh::DofMap::is_constrained_dof ( const dof_id_type  dof) const
inline
Returns
true if the degree of freedom dof is constrained, false otherwise.

Definition at line 1432 of file dof_map.h.

References _dof_constraints.

Referenced by libMesh::FEGenericBase< T >::compute_periodic_constraints(), libMesh::FEGenericBase< T >::compute_proj_constraints(), extract_local_vector(), and find_connected_dofs().

1433 {
1434  if (_dof_constraints.count(dof))
1435  return true;
1436 
1437  return false;
1438 }
bool libMesh::DofMap::is_constrained_node ( const Node node) const
inline
Returns
true if the Node is constrained, false otherwise.

Definition at line 1415 of file dof_map.h.

References _node_constraints.

1420 {
1421 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1422  if (_node_constraints.count(node))
1423  return true;
1424 #endif
1425 
1426  return false;
1427 }
bool libMesh::DofMap::is_periodic_boundary ( const boundary_id_type  boundaryid) const
Returns
true if the boundary given by boundaryid is periodic, false otherwise

Definition at line 193 of file dof_map.C.

References _periodic_boundaries.

194 {
195  if (_periodic_boundaries->count(boundaryid) != 0)
196  return true;
197 
198  return false;
199 }
dof_id_type libMesh::DofMap::last_dof ( const processor_id_type  proc) const
inline

Returns the last dof index that is local to processor 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 475 of file dof_map.h.

References _end_df.

475  {
476  libmesh_deprecated();
477  libmesh_assert_less (proc, _end_df.size());
478  return libmesh_cast_int<dof_id_type>(_end_df[proc] - 1);
479  }
dof_id_type libMesh::DofMap::last_dof ( ) const
inline

Definition at line 481 of file dof_map.h.

References libMesh::ParallelObject::processor_id().

482  { return this->last_dof(this->processor_id()); }
std::pair< Real, Real > libMesh::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 2051 of file dof_map_constraints.C.

References std::abs(), libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::NumericVector< T >::closed(), libMesh::NumericVector< T >::first_local_index(), libMesh::for(), libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::NumericVector< T >::last_local_index(), libMesh::libmesh_assert(), libMesh::DenseMatrixBase< T >::m(), std::max(), mesh, libMesh::DenseMatrixBase< T >::n(), libMesh::Real, and libMesh::System::solution.

2053 {
2054  if (!v)
2055  v = system.solution.get();
2056  NumericVector<Number> &vec = *v;
2057 
2058  // We'll assume the vector is closed
2059  libmesh_assert (vec.closed());
2060 
2061  Real max_absolute_error = 0., max_relative_error = 0.;
2062 
2063  const MeshBase &mesh = system.get_mesh();
2064 
2065  libmesh_assert_equal_to (this, &(system.get_dof_map()));
2066 
2067  // indices on each element
2068  std::vector<dof_id_type> local_dof_indices;
2069 
2072  const MeshBase::const_element_iterator elem_end =
2074 
2075  for ( ; elem_it != elem_end; ++elem_it)
2076  {
2077  const Elem* elem = *elem_it;
2078 
2079  this->dof_indices(elem, local_dof_indices);
2080  std::vector<dof_id_type> raw_dof_indices = local_dof_indices;
2081 
2082  // Constraint matrix for each element
2084 
2085  this->build_constraint_matrix (C, local_dof_indices);
2086 
2087  // Continue if the element is unconstrained
2088  if (!C.m())
2089  continue;
2090 
2091  libmesh_assert_equal_to (C.m(), raw_dof_indices.size());
2092  libmesh_assert_equal_to (C.n(), local_dof_indices.size());
2093 
2094  for (unsigned int i=0; i!=C.m(); ++i)
2095  {
2096  // Recalculate any constrained dof owned by this processor
2097  dof_id_type global_dof = raw_dof_indices[i];
2098  if (this->is_constrained_dof(global_dof) &&
2099  global_dof >= vec.first_local_index() &&
2100  global_dof < vec.last_local_index())
2101  {
2102  DofConstraints::const_iterator
2103  pos = _dof_constraints.find(global_dof);
2104 
2105  libmesh_assert (pos != _dof_constraints.end());
2106 
2107  Number exact_value = 0;
2108  DofConstraintValueMap::const_iterator rhsit =
2109  _primal_constraint_values.find(global_dof);
2110  if (rhsit != _primal_constraint_values.end())
2111  exact_value = rhsit->second;
2112 
2113  for (unsigned int j=0; j!=C.n(); ++j)
2114  {
2115  if (local_dof_indices[j] != global_dof)
2116  exact_value += C(i,j) *
2117  vec(local_dof_indices[j]);
2118  }
2119 
2120  max_absolute_error = std::max(max_absolute_error,
2121  std::abs(vec(global_dof) - exact_value));
2122  max_relative_error = std::max(max_relative_error,
2123  std::abs(vec(global_dof) - exact_value)
2124  / std::abs(exact_value));
2125  }
2126  }
2127  }
2128 
2129  return std::pair<Real, Real>(max_absolute_error, max_relative_error);
2130 }
dof_id_type libMesh::DofMap::n_constrained_dofs ( ) const
Returns
the total number of constrained degrees of freedom in the problem.

Definition at line 919 of file dof_map_constraints.C.

References libMesh::comm.

920 {
921  parallel_object_only();
922 
923  dof_id_type nc_dofs = this->n_local_constrained_dofs();
924  this->comm().sum(nc_dofs);
925  return nc_dofs;
926 }
dof_id_type libMesh::DofMap::n_constrained_nodes ( ) const
inline
Returns
the total number of constrained Nodes in the mesh.

Definition at line 578 of file dof_map.h.

References _node_constraints.

579  { return libmesh_cast_int<dof_id_type>(_node_constraints.size()); }
dof_id_type libMesh::DofMap::n_dofs ( ) const
inline
Returns
the total number of degrees of freedom in the problem.

Definition at line 428 of file dof_map.h.

References _n_dfs.

Referenced by libMesh::SparsityPattern::Build::join(), libMesh::UnsteadySolver::old_nonlinear_solution(), libMesh::SparsityPattern::Build::parallel_sync(), and SCALAR_dof_indices().

428 { return _n_dfs; }
dof_id_type libMesh::DofMap::n_dofs_on_processor ( const processor_id_type  proc) const
inline

Returns the number of degrees of freedom on partition proc.

Definition at line 444 of file dof_map.h.

References _end_df, and _first_df.

Referenced by libMesh::SparsityPattern::Build::join(), n_local_dofs(), libMesh::SparsityPattern::Build::operator()(), and libMesh::SparsityPattern::Build::parallel_sync().

444  {
445  libmesh_assert_less (proc, _first_df.size());
446  return libmesh_cast_int<dof_id_type>(_end_df[proc] - _first_df[proc]);
447  }
dof_id_type libMesh::DofMap::n_local_constrained_dofs ( ) const
Returns
the number of constrained degrees of freedom on this processor.

Definition at line 929 of file dof_map_constraints.C.

930 {
931  const DofConstraints::const_iterator lower =
932  _dof_constraints.lower_bound(this->first_dof()),
933  upper =
934  _dof_constraints.upper_bound(this->end_dof()-1);
935 
936  return std::distance(lower, upper);
937 }
dof_id_type libMesh::DofMap::n_local_dofs ( ) const
inline
Returns
the number of degrees of freedom on this processor.

Definition at line 438 of file dof_map.h.

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

439  { return this->n_dofs_on_processor (this->processor_id()); }
static unsigned int libMesh::ReferenceCounter::n_objects ( )
inlinestaticinherited

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

Definition at line 79 of file reference_counter.h.

References libMesh::ReferenceCounter::_n_objects.

80  { return _n_objects; }
dof_id_type libMesh::DofMap::n_old_dofs ( ) const
inline
Returns
the total number of degrees of freedom on old_dof_objects

Definition at line 991 of file dof_map.h.

References _n_old_dfs.

Referenced by SCALAR_dof_indices().

991 { return _n_old_dfs; }
processor_id_type libMesh::ParallelObject::n_processors ( ) const
inlineinherited
Returns
the number of processors in the group.

Definition at line 92 of file parallel_object.h.

References libMesh::ParallelObject::_communicator, and libMesh::Parallel::Communicator::size().

Referenced by libMesh::ParmetisPartitioner::_do_repartition(), libMesh::ParallelMesh::add_elem(), libMesh::ParallelMesh::add_node(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::ParmetisPartitioner::assign_partitioning(), libMesh::ParallelMesh::assign_unique_ids(), libMesh::AztecLinearSolver< T >::AztecLinearSolver(), libMesh::ParallelMesh::clear(), libMesh::Nemesis_IO_Helper::compute_border_node_ids(), libMesh::Nemesis_IO_Helper::construct_nemesis_filename(), libMesh::UnstructuredMesh::create_pid_mesh(), distribute_dofs(), distribute_local_dofs_node_major(), distribute_local_dofs_var_major(), libMesh::EnsightIO::EnsightIO(), libMesh::MeshBase::get_info(), libMesh::EquationSystems::init(), libMesh::SystemSubsetBySubdomain::init(), libMesh::ParmetisPartitioner::initialize(), libMesh::Nemesis_IO_Helper::initialize(), libMesh::MeshTools::libmesh_assert_valid_dof_ids(), libMesh::MeshTools::libmesh_assert_valid_procids< Elem >(), libMesh::MeshTools::libmesh_assert_valid_procids< Node >(), libMesh::MeshTools::libmesh_assert_valid_refinement_flags(), libMesh::MeshBase::n_active_elem_on_proc(), libMesh::MeshBase::n_elem_on_proc(), libMesh::MeshBase::n_nodes_on_proc(), libMesh::Partitioner::partition(), libMesh::MeshBase::partition(), libMesh::Partitioner::partition_unpartitioned_elements(), libMesh::PetscLinearSolver< T >::PetscLinearSolver(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::MeshTools::processor_bounding_box(), libMesh::System::project_vector(), libMesh::Nemesis_IO::read(), libMesh::CheckpointIO::read(), libMesh::UnstructuredMesh::read(), libMesh::System::read_parallel_data(), libMesh::System::read_SCALAR_dofs(), libMesh::System::read_serialized_blocked_dof_objects(), libMesh::System::read_serialized_vector(), libMesh::Partitioner::repartition(), libMesh::Partitioner::set_node_processor_ids(), set_nonlocal_dof_objects(), libMesh::BoundaryInfo::sync(), libMesh::ParallelMesh::update_parallel_id_counts(), libMesh::CheckpointIO::write(), libMesh::GMVIO::write_binary(), libMesh::GMVIO::write_discontinuous_gmv(), libMesh::System::write_parallel_data(), libMesh::System::write_SCALAR_dofs(), libMesh::XdrIO::write_serialized_bcs(), libMesh::System::write_serialized_blocked_dof_objects(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), and libMesh::XdrIO::write_serialized_nodesets().

93  { return libmesh_cast_int<processor_id_type>(_communicator.size()); }
dof_id_type libMesh::DofMap::n_SCALAR_dofs ( ) const
inline
Returns
the number of SCALAR dofs.

Definition at line 433 of file dof_map.h.

References _n_SCALAR_dofs.

Referenced by SCALAR_dof_indices().

433 { return _n_SCALAR_dofs; }
unsigned int libMesh::DofMap::n_variable_groups ( ) 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 387 of file dof_map.h.

References _variable_groups.

Referenced by distribute_local_dofs_node_major(), distribute_local_dofs_var_major(), has_blocked_representation(), reinit(), and set_nonlocal_dof_objects().

388  { return libmesh_cast_int<unsigned int>(_variable_groups.size()); }
unsigned int libMesh::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 395 of file dof_map.h.

References _variables.

Referenced by block_size(), DMLibMeshSetSystem(), dof_indices(), has_blocked_representation(), old_dof_indices(), libMesh::SparsityPattern::Build::operator()(), SCALAR_dof_indices(), and use_coupled_neighbor_dofs().

396  { return libmesh_cast_int<unsigned int>(_variables.size()); }
NodeConstraints::const_iterator libMesh::DofMap::node_constraint_rows_begin ( ) const
inline

Returns an iterator pointing to the first Node constraint row

Definition at line 659 of file dof_map.h.

References _node_constraints.

660  { return _node_constraints.begin(); }
NodeConstraints::const_iterator libMesh::DofMap::node_constraint_rows_end ( ) const
inline

Returns an iterator pointing just past the last Node constraint row

Definition at line 665 of file dof_map.h.

References _node_constraints.

666  { return _node_constraints.end(); }
DofObject * libMesh::DofMap::node_ptr ( MeshBase mesh,
dof_id_type  i 
) const
private

An adapter function that returns Node pointers by index

Definition at line 269 of file dof_map.C.

References libMesh::MeshBase::node_ptr().

Referenced by distribute_dofs().

270 {
271  return mesh.node_ptr(i);
272 }
void libMesh::DofMap::old_dof_indices ( const Elem *const  elem,
std::vector< dof_id_type > &  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 1886 of file dof_map.C.

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

Referenced by libMesh::FEGenericBase< T >::coarsened_dof_values().

1889 {
1890  START_LOG("old_dof_indices()", "DofMap");
1891 
1892  libmesh_assert(elem);
1893  libmesh_assert(elem->old_dof_object);
1894 
1895 
1896  const unsigned int n_nodes = elem->n_nodes();
1897  const ElemType type = elem->type();
1898  const unsigned int sys_num = this->sys_number();
1899  const unsigned int n_vars = this->n_variables();
1900  const unsigned int dim = elem->dim();
1901 
1902  // Clear the DOF indices vector.
1903  di.clear();
1904 
1905 #ifdef DEBUG
1906  // Check that sizes match
1907  unsigned int tot_size = 0;
1908 #endif
1909 
1910  // Create a vector to indicate which
1911  // SCALAR variables have been requested
1912  std::vector<unsigned int> SCALAR_var_numbers;
1913  SCALAR_var_numbers.clear();
1914 
1915  // Get the dof numbers
1916  for (unsigned int v=0; v<n_vars; v++)
1917  if ((v == vn) || (vn == libMesh::invalid_uint))
1918  {
1919  if(this->variable(v).type().family == SCALAR)
1920  {
1921  // We asked for this variable, so add it to the vector.
1922  SCALAR_var_numbers.push_back(v);
1923 
1924 #ifdef DEBUG
1925  FEType fe_type = this->variable_type(v);
1926  tot_size += FEInterface::n_dofs(dim,
1927  fe_type,
1928  type);
1929 #endif
1930  }
1931  else
1932  if (this->variable(v).active_on_subdomain(elem->subdomain_id()))
1933  { // Do this for all the variables if one was not specified
1934  // or just for the specified variable
1935 
1936  // Increase the polynomial order on p refined elements,
1937  // but make sure you get the right polynomial order for
1938  // the OLD degrees of freedom
1939  int p_adjustment = 0;
1940  if (elem->p_refinement_flag() == Elem::JUST_REFINED)
1941  {
1942  libmesh_assert_greater (elem->p_level(), 0);
1943  p_adjustment = -1;
1944  }
1945  else if (elem->p_refinement_flag() == Elem::JUST_COARSENED)
1946  {
1947  p_adjustment = 1;
1948  }
1949  FEType fe_type = this->variable_type(v);
1950  fe_type.order = static_cast<Order>(fe_type.order +
1951  elem->p_level() +
1952  p_adjustment);
1953 
1954  const bool extra_hanging_dofs =
1956 
1957 #ifdef DEBUG
1958  tot_size += FEInterface::n_dofs(dim,
1959  fe_type,
1960  type);
1961 #endif
1962 
1963  // Get the node-based DOF numbers
1964  for (unsigned int n=0; n<n_nodes; n++)
1965  {
1966  const Node* node = elem->get_node(n);
1967 
1968  // There is a potential problem with h refinement. Imagine a
1969  // quad9 that has a linear FE on it. Then, on the hanging side,
1970  // it can falsely identify a DOF at the mid-edge node. This is why
1971  // we call FEInterface instead of node->n_comp() directly.
1972  const unsigned int nc = FEInterface::n_dofs_at_node (dim,
1973  fe_type,
1974  type,
1975  n);
1976  libmesh_assert(node->old_dof_object);
1977 
1978  // If this is a non-vertex on a hanging node with extra
1979  // degrees of freedom, we use the non-vertex dofs (which
1980  // come in reverse order starting from the end, to
1981  // simplify p refinement)
1982  if (extra_hanging_dofs && !elem->is_vertex(n))
1983  {
1984  const int dof_offset =
1985  node->old_dof_object->n_comp(sys_num,v) - nc;
1986 
1987  // We should never have fewer dofs than necessary on a
1988  // node unless we're getting indices on a parent element
1989  // or a just-coarsened element
1990  if (dof_offset < 0)
1991  {
1992  libmesh_assert(!elem->active() || elem->refinement_flag() ==
1994  di.resize(di.size() + nc, DofObject::invalid_id);
1995  }
1996  else
1997  for (int i=node->old_dof_object->n_comp(sys_num,v)-1;
1998  i>=dof_offset; i--)
1999  {
2000  libmesh_assert_not_equal_to (node->old_dof_object->dof_number(sys_num,v,i),
2002  di.push_back(node->old_dof_object->dof_number(sys_num,v,i));
2003  }
2004  }
2005  // If this is a vertex or an element without extra hanging
2006  // dofs, our dofs come in forward order coming from the
2007  // beginning
2008  else
2009  for (unsigned int i=0; i<nc; i++)
2010  {
2011  libmesh_assert_not_equal_to (node->old_dof_object->dof_number(sys_num,v,i),
2013  di.push_back(node->old_dof_object->dof_number(sys_num,v,i));
2014  }
2015  }
2016 
2017  // If there are any element-based DOF numbers, get them
2018  const unsigned int nc = FEInterface::n_dofs_per_elem(dim,
2019  fe_type,
2020  type);
2021 
2022  // We should never have fewer dofs than necessary on an
2023  // element unless we're getting indices on a parent element
2024  // or a just-coarsened element
2025  if (nc != 0)
2026  {
2027  if (elem->old_dof_object->n_systems() > sys_num &&
2028  nc <= elem->old_dof_object->n_comp(sys_num,v))
2029  {
2030  libmesh_assert(elem->old_dof_object);
2031 
2032  for (unsigned int i=0; i<nc; i++)
2033  {
2034  libmesh_assert_not_equal_to (elem->old_dof_object->dof_number(sys_num,v,i),
2036 
2037  di.push_back(elem->old_dof_object->dof_number(sys_num,v,i));
2038  }
2039  }
2040  else
2041  {
2042  libmesh_assert(!elem->active() || fe_type.family == LAGRANGE ||
2043  elem->refinement_flag() == Elem::JUST_COARSENED);
2044  di.resize(di.size() + nc, DofObject::invalid_id);
2045  }
2046  }
2047  }
2048  } // end loop over variables
2049 
2050  // Finally append any SCALAR dofs that we asked for.
2051  std::vector<dof_id_type> di_new;
2052  std::vector<unsigned int>::iterator it = SCALAR_var_numbers.begin();
2053  std::vector<unsigned int>::const_iterator it_end = SCALAR_var_numbers.end();
2054  for( ; it != it_end; ++it)
2055  {
2056  this->SCALAR_dof_indices(di_new,*it,true);
2057  di.insert( di.end(), di_new.begin(), di_new.end());
2058  }
2059 
2060 #ifdef DEBUG
2061  libmesh_assert_equal_to (tot_size, di.size());
2062 #endif
2063 
2064  STOP_LOG("old_dof_indices()", "DofMap");
2065 }
void libMesh::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. Also calls any user-provided methods for adding to the send list.

Definition at line 1428 of file dof_map.C.

References _augment_send_list, _extra_send_list_context, _extra_send_list_function, _send_list, libMesh::DofMap::AugmentSendList::augment_send_list(), libMesh::out, libMesh::START_LOG(), libMesh::STOP_LOG(), and swap().

Referenced by libMesh::EquationSystems::allgather(), and libMesh::EquationSystems::reinit().

1429 {
1430  START_LOG("prepare_send_list()", "DofMap");
1431 
1432  // Check to see if we have any extra stuff to add to the send_list
1434  {
1435  if (_augment_send_list)
1436  {
1437  libmesh_here();
1438  libMesh::out << "WARNING: You have specified both an extra send list function and object.\n"
1439  << " Are you sure this is what you meant to do??"
1440  << std::endl;
1441  }
1442 
1444  }
1445 
1446  if (_augment_send_list)
1448 
1449  // First sort the send list. After this
1450  // duplicated elements will be adjacent in the
1451  // vector
1452  std::sort(_send_list.begin(), _send_list.end());
1453 
1454  // Now use std::unique to remove duplicate entries
1455  std::vector<dof_id_type>::iterator new_end =
1456  std::unique (_send_list.begin(), _send_list.end());
1457 
1458  // Remove the end of the send_list. Use the "swap trick"
1459  // from Effective STL
1460  std::vector<dof_id_type> (_send_list.begin(), new_end).swap (_send_list);
1461 
1462  STOP_LOG("prepare_send_list()", "DofMap");
1463 }
void libMesh::DofMap::print_dof_constraints ( std::ostream &  os = libMesh::out,
bool  print_nonlocal = false 
) const

Prints (from processor 0) all DoF and Node constraints. If print_nonlocal is true, then each constraint is printed once for each processor that knows about it, which may be useful for ParallelMesh debugging.

Definition at line 1164 of file dof_map_constraints.C.

References libMesh::comm, libMesh::n_processors(), and libMesh::processor_id().

1166 {
1167  parallel_object_only();
1168 
1169  std::string local_constraints =
1170  this->get_local_constraints(print_nonlocal);
1171 
1172  if (this->processor_id())
1173  {
1174  this->comm().send(0, local_constraints);
1175  }
1176  else
1177  {
1178  os << "Processor 0:\n";
1179  os << local_constraints;
1180 
1181  for (processor_id_type i=1; i<this->n_processors(); ++i)
1182  {
1183  this->comm().receive(i, local_constraints);
1184  os << "Processor " << i << ":\n";
1185  os << local_constraints;
1186  }
1187  }
1188 }
void libMesh::ReferenceCounter::print_info ( std::ostream &  out = libMesh::out)
staticinherited

Prints the reference information, by default to libMesh::out.

Definition at line 88 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter, and libMesh::ReferenceCounter::get_info().

89 {
91 }
void libMesh::DofMap::print_info ( std::ostream &  os = libMesh::out) const

Prints summary info about the sparsity bandwidth and constraints.

Definition at line 2824 of file dof_map.C.

2825 {
2826  os << this->get_info();
2827 }
void libMesh::DofMap::process_constraints ( MeshBase mesh)

Postprocesses any constrained degrees of freedom 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 2927 of file dof_map_constraints.C.

References libMesh::for(), libMesh::MeshBase::is_serial(), libMesh::libmesh_assert(), and libMesh::Real.

Referenced by libMesh::EquationSystems::allgather(), and libMesh::EquationSystems::reinit().

2928 {
2929  // With a parallelized Mesh, we've computed our local constraints,
2930  // but they may depend on non-local constraints that we'll need to
2931  // take into account.
2932  if (!mesh.is_serial())
2933  this->allgather_recursive_constraints(mesh);
2934 
2935  // Create a set containing the DOFs we already depend on
2936  typedef std::set<dof_id_type> RCSet;
2937  RCSet unexpanded_set;
2938 
2939  for (DofConstraints::iterator i = _dof_constraints.begin();
2940  i != _dof_constraints.end(); ++i)
2941  unexpanded_set.insert(i->first);
2942 
2943  while (!unexpanded_set.empty())
2944  for (RCSet::iterator i = unexpanded_set.begin();
2945  i != unexpanded_set.end(); /* nothing */)
2946  {
2947  // If the DOF is constrained
2948  DofConstraints::iterator
2949  pos = _dof_constraints.find(*i);
2950 
2951  libmesh_assert (pos != _dof_constraints.end());
2952 
2953  DofConstraintRow& constraint_row = pos->second;
2954 
2955  DofConstraintValueMap::iterator rhsit =
2956  _primal_constraint_values.find(*i);
2957  Number constraint_rhs = (rhsit == _primal_constraint_values.end()) ?
2958  0 : rhsit->second;
2959 
2960  std::vector<dof_id_type> constraints_to_expand;
2961 
2962  for (DofConstraintRow::const_iterator
2963  it=constraint_row.begin(); it != constraint_row.end();
2964  ++it)
2965  if (it->first != *i && this->is_constrained_dof(it->first))
2966  {
2967  unexpanded_set.insert(it->first);
2968  constraints_to_expand.push_back(it->first);
2969  }
2970 
2971  for (std::size_t j = 0; j != constraints_to_expand.size();
2972  ++j)
2973  {
2974  dof_id_type expandable = constraints_to_expand[j];
2975 
2976  const Real this_coef = constraint_row[expandable];
2977 
2978  DofConstraints::const_iterator
2979  subpos = _dof_constraints.find(expandable);
2980 
2981  libmesh_assert (subpos != _dof_constraints.end());
2982 
2983  const DofConstraintRow& subconstraint_row = subpos->second;
2984 
2985  for (DofConstraintRow::const_iterator
2986  it=subconstraint_row.begin();
2987  it != subconstraint_row.end(); ++it)
2988  {
2989  constraint_row[it->first] += it->second * this_coef;
2990  }
2991  DofConstraintValueMap::const_iterator subrhsit =
2992  _primal_constraint_values.find(expandable);
2993  if (subrhsit != _primal_constraint_values.end())
2994  constraint_rhs += subrhsit->second * this_coef;
2995 
2996  constraint_row.erase(expandable);
2997  }
2998 
2999  if (rhsit == _primal_constraint_values.end())
3000  {
3001  if (constraint_rhs != Number(0))
3002  _primal_constraint_values[*i] = constraint_rhs;
3003  else
3004  _primal_constraint_values.erase(*i);
3005  }
3006  else
3007  {
3008  if (constraint_rhs != Number(0))
3009  rhsit->second = constraint_rhs;
3010  else
3011  _primal_constraint_values.erase(rhsit);
3012  }
3013 
3014  if (constraints_to_expand.empty())
3015  unexpanded_set.erase(i++);
3016  else
3017  i++;
3018  }
3019 
3020  // In parallel we can't guarantee that nodes/dofs which constrain
3021  // others are on processors which are aware of that constraint, yet
3022  // we need such awareness for sparsity pattern generation. So send
3023  // other processors any constraints they might need to know about.
3024  if (!mesh.is_serial())
3025  this->scatter_constraints(mesh);
3026 
3027  // Now that we have our root constraint dependencies sorted out, add
3028  // them to the send_list
3030 }
processor_id_type libMesh::ParallelObject::processor_id ( ) const
inlineinherited
Returns
the rank of this processor in the group.

Definition at line 98 of file parallel_object.h.

References libMesh::ParallelObject::_communicator, and libMesh::Parallel::Communicator::rank().

Referenced by libMesh::MetisPartitioner::_do_partition(), libMesh::EquationSystems::_read_impl(), libMesh::SerialMesh::active_local_elements_begin(), libMesh::ParallelMesh::active_local_elements_begin(), libMesh::SerialMesh::active_local_elements_end(), libMesh::ParallelMesh::active_local_elements_end(), libMesh::SerialMesh::active_local_subdomain_elements_begin(), libMesh::ParallelMesh::active_local_subdomain_elements_begin(), libMesh::SerialMesh::active_local_subdomain_elements_end(), libMesh::ParallelMesh::active_local_subdomain_elements_end(), libMesh::SerialMesh::active_not_local_elements_begin(), libMesh::ParallelMesh::active_not_local_elements_begin(), libMesh::SerialMesh::active_not_local_elements_end(), libMesh::ParallelMesh::active_not_local_elements_end(), libMesh::ParallelMesh::add_elem(), add_neighbors_to_send_list(), libMesh::ParallelMesh::add_node(), libMesh::UnstructuredMesh::all_second_order(), libMesh::ParmetisPartitioner::assign_partitioning(), libMesh::EquationSystems::build_discontinuous_solution_vector(), libMesh::Nemesis_IO_Helper::build_element_and_node_maps(), libMesh::ParmetisPartitioner::build_graph(), libMesh::InfElemBuilder::build_inf_elem(), build_sparsity(), libMesh::ParallelMesh::clear(), libMesh::ExodusII_IO_Helper::close(), libMesh::Nemesis_IO_Helper::compute_border_node_ids(), libMesh::Nemesis_IO_Helper::compute_communication_map_parameters(), libMesh::Nemesis_IO_Helper::compute_internal_and_border_elems_and_internal_nodes(), libMesh::Nemesis_IO_Helper::compute_node_communication_maps(), libMesh::Nemesis_IO_Helper::compute_num_global_elem_blocks(), libMesh::Nemesis_IO_Helper::compute_num_global_nodesets(), libMesh::Nemesis_IO_Helper::compute_num_global_sidesets(), libMesh::Nemesis_IO_Helper::construct_nemesis_filename(), libMesh::ExodusII_IO_Helper::create(), distribute_dofs(), distribute_local_dofs_node_major(), distribute_local_dofs_var_major(), end_dof(), end_old_dof(), libMesh::EnsightIO::EnsightIO(), libMesh::UnstructuredMesh::find_neighbors(), first_dof(), first_old_dof(), libMesh::Nemesis_IO_Helper::get_cmap_params(), libMesh::Nemesis_IO_Helper::get_eb_info_global(), libMesh::Nemesis_IO_Helper::get_elem_cmap(), libMesh::Nemesis_IO_Helper::get_elem_map(), libMesh::MeshBase::get_info(), libMesh::Nemesis_IO_Helper::get_init_global(), libMesh::Nemesis_IO_Helper::get_init_info(), libMesh::Nemesis_IO_Helper::get_loadbal_param(), libMesh::Nemesis_IO_Helper::get_node_cmap(), libMesh::Nemesis_IO_Helper::get_node_map(), libMesh::Nemesis_IO_Helper::get_ns_param_global(), libMesh::Nemesis_IO_Helper::get_ss_param_global(), libMesh::MeshFunction::gradient(), libMesh::MeshFunction::hessian(), libMesh::SystemSubsetBySubdomain::init(), libMesh::ParmetisPartitioner::initialize(), libMesh::ExodusII_IO_Helper::initialize(), libMesh::ExodusII_IO_Helper::initialize_discontinuous(), libMesh::ExodusII_IO_Helper::initialize_element_variables(), libMesh::ExodusII_IO_Helper::initialize_global_variables(), libMesh::ExodusII_IO_Helper::initialize_nodal_variables(), libMesh::SparsityPattern::Build::join(), last_dof(), libMesh::MeshTools::libmesh_assert_valid_procids< Elem >(), libMesh::MeshTools::libmesh_assert_valid_procids< Node >(), libMesh::SerialMesh::local_elements_begin(), libMesh::ParallelMesh::local_elements_begin(), libMesh::SerialMesh::local_elements_end(), libMesh::ParallelMesh::local_elements_end(), libMesh::SerialMesh::local_level_elements_begin(), libMesh::ParallelMesh::local_level_elements_begin(), libMesh::SerialMesh::local_level_elements_end(), libMesh::ParallelMesh::local_level_elements_end(), libMesh::SerialMesh::local_nodes_begin(), libMesh::ParallelMesh::local_nodes_begin(), libMesh::SerialMesh::local_nodes_end(), libMesh::ParallelMesh::local_nodes_end(), libMesh::SerialMesh::local_not_level_elements_begin(), libMesh::ParallelMesh::local_not_level_elements_begin(), libMesh::SerialMesh::local_not_level_elements_end(), libMesh::ParallelMesh::local_not_level_elements_end(), libMesh::MeshRefinement::make_coarsening_compatible(), libMesh::MeshBase::n_active_local_elem(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::BoundaryInfo::n_edge_conds(), n_local_dofs(), libMesh::System::n_local_dofs(), libMesh::MeshBase::n_local_elem(), libMesh::MeshBase::n_local_nodes(), libMesh::BoundaryInfo::n_nodeset_conds(), libMesh::SerialMesh::not_local_elements_begin(), libMesh::ParallelMesh::not_local_elements_begin(), libMesh::SerialMesh::not_local_elements_end(), libMesh::ParallelMesh::not_local_elements_end(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::SparsityPattern::Build::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::MeshFunction::operator()(), libMesh::ParallelMesh::ParallelMesh(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::System::project_vector(), libMesh::Nemesis_IO_Helper::put_cmap_params(), libMesh::Nemesis_IO_Helper::put_elem_cmap(), libMesh::Nemesis_IO_Helper::put_elem_map(), libMesh::Nemesis_IO_Helper::put_loadbal_param(), libMesh::Nemesis_IO_Helper::put_node_cmap(), libMesh::Nemesis_IO_Helper::put_node_map(), libMesh::Nemesis_IO::read(), libMesh::CheckpointIO::read(), libMesh::XdrIO::read(), libMesh::UnstructuredMesh::read(), libMesh::CheckpointIO::read_connectivity(), libMesh::ExodusII_IO_Helper::read_elem_num_map(), libMesh::System::read_header(), libMesh::System::read_legacy_data(), libMesh::ExodusII_IO_Helper::read_node_num_map(), libMesh::System::read_parallel_data(), libMesh::System::read_SCALAR_dofs(), libMesh::XdrIO::read_serialized_bc_names(), libMesh::XdrIO::read_serialized_bcs(), libMesh::System::read_serialized_blocked_dof_objects(), libMesh::XdrIO::read_serialized_connectivity(), libMesh::System::read_serialized_data(), libMesh::XdrIO::read_serialized_nodes(), libMesh::XdrIO::read_serialized_nodesets(), libMesh::XdrIO::read_serialized_subdomain_names(), libMesh::System::read_serialized_vector(), libMesh::System::read_serialized_vectors(), libMesh::MeshData::read_xdr(), libMesh::Partitioner::set_node_processor_ids(), set_nonlocal_dof_objects(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::BoundaryInfo::sync(), libMesh::MeshTools::total_weight(), libMesh::ParallelMesh::update_parallel_id_counts(), libMesh::MeshTools::weight(), libMesh::ExodusII_IO::write(), libMesh::CheckpointIO::write(), libMesh::XdrIO::write(), libMesh::UnstructuredMesh::write(), libMesh::EquationSystems::write(), libMesh::GMVIO::write_discontinuous_gmv(), libMesh::ExodusII_IO::write_element_data(), libMesh::ExodusII_IO_Helper::write_element_values(), libMesh::ExodusII_IO_Helper::write_elements(), libMesh::ExodusII_IO_Helper::write_elements_discontinuous(), libMesh::ExodusII_IO::write_global_data(), libMesh::ExodusII_IO_Helper::write_global_values(), libMesh::System::write_header(), libMesh::ExodusII_IO::write_information_records(), libMesh::ExodusII_IO_Helper::write_information_records(), libMesh::ExodusII_IO_Helper::write_nodal_coordinates(), libMesh::ExodusII_IO_Helper::write_nodal_coordinates_discontinuous(), libMesh::UCDIO::write_nodal_data(), libMesh::ExodusII_IO::write_nodal_data(), libMesh::ExodusII_IO::write_nodal_data_discontinuous(), libMesh::ExodusII_IO_Helper::write_nodal_values(), libMesh::ExodusII_IO_Helper::write_nodesets(), libMesh::Nemesis_IO_Helper::write_nodesets(), libMesh::System::write_parallel_data(), libMesh::System::write_SCALAR_dofs(), libMesh::XdrIO::write_serialized_bc_names(), libMesh::XdrIO::write_serialized_bcs(), libMesh::System::write_serialized_blocked_dof_objects(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::System::write_serialized_data(), libMesh::XdrIO::write_serialized_nodes(), libMesh::XdrIO::write_serialized_nodesets(), libMesh::XdrIO::write_serialized_subdomain_names(), libMesh::System::write_serialized_vector(), libMesh::System::write_serialized_vectors(), libMesh::ExodusII_IO_Helper::write_sidesets(), libMesh::Nemesis_IO_Helper::write_sidesets(), libMesh::ExodusII_IO::write_timestep(), and libMesh::ExodusII_IO_Helper::write_timestep().

99  { return libmesh_cast_int<processor_id_type>(_communicator.rank()); }
void libMesh::DofMap::reinit ( MeshBase mesh)

Reinitialize the underlying data strucures conformal to the current mesh.

Definition at line 432 of file dof_map.C.

References _n_SCALAR_dofs, libMesh::MeshBase::active_elements_begin(), libMesh::MeshBase::active_elements_end(), libMesh::Variable::active_on_subdomain(), libMesh::dim, libMesh::Elem::dim(), libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), libMesh::Utility::enum_to_string(), libMesh::err, libMesh::FEInterface::extra_hanging_dofs(), libMesh::FEType::family, libMesh::Elem::get_node(), libMesh::DofObject::has_dofs(), invalidate_dofs(), libMesh::MeshBase::is_prepared(), libMesh::Elem::is_vertex(), libMesh::Elem::JUST_REFINED, libMesh::libmesh_assert(), std::max(), libMesh::FEInterface::max_order(), libMesh::DofObject::n_comp_group(), libMesh::FEInterface::n_dofs_at_node(), libMesh::FEInterface::n_dofs_per_elem(), libMesh::Elem::n_nodes(), n_variable_groups(), libMesh::VariableGroup::n_variables(), libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), libMesh::DofObject::old_dof_object, libMesh::FEType::order, libMesh::Elem::p_level(), libMesh::Elem::refinement_flag(), libMeshEnums::SCALAR, libMesh::DofObject::set_n_comp_group(), libMesh::DofObject::set_old_dof_object(), libMesh::Elem::set_p_level(), libMesh::DofObject::set_vg_dof_base(), libMesh::START_LOG(), libMesh::STOP_LOG(), libMesh::Elem::subdomain_id(), sys_number(), libMesh::Variable::type(), libMesh::Elem::type(), variable_group(), and libMesh::DofObject::vg_dof_base().

Referenced by distribute_dofs().

433 {
434  libmesh_assert (mesh.is_prepared());
435 
436  START_LOG("reinit()", "DofMap");
437 
438  const unsigned int
439  sys_num = this->sys_number(),
440  n_var_groups = this->n_variable_groups();
441 
442  // The DofObjects need to know how many variable groups we have, and
443  // how many variables there are in each group.
444  std::vector<unsigned int> n_vars_per_group; n_vars_per_group.reserve (n_var_groups);
445 
446  for (unsigned int vg=0; vg<n_var_groups; vg++)
447  n_vars_per_group.push_back (this->variable_group(vg).n_variables());
448 
449 #ifdef LIBMESH_ENABLE_AMR
450 
451  //------------------------------------------------------------
452  // Clear the old_dof_objects for all the nodes
453  // and elements so that we can overwrite them
454  {
455  MeshBase::node_iterator node_it = mesh.nodes_begin();
456  const MeshBase::node_iterator node_end = mesh.nodes_end();
457 
458  for ( ; node_it != node_end; ++node_it)
459  {
460  (*node_it)->clear_old_dof_object();
461  libmesh_assert (!(*node_it)->old_dof_object);
462  }
463 
464  MeshBase::element_iterator elem_it = mesh.elements_begin();
465  const MeshBase::element_iterator elem_end = mesh.elements_end();
466 
467  for ( ; elem_it != elem_end; ++elem_it)
468  {
469  (*elem_it)->clear_old_dof_object();
470  libmesh_assert (!(*elem_it)->old_dof_object);
471  }
472  }
473 
474 
475  //------------------------------------------------------------
476  // Set the old_dof_objects for the elements that
477  // weren't just created, if these old dof objects
478  // had variables
479  {
480  MeshBase::element_iterator elem_it = mesh.elements_begin();
481  const MeshBase::element_iterator elem_end = mesh.elements_end();
482 
483  for ( ; elem_it != elem_end; ++elem_it)
484  {
485  Elem* elem = *elem_it;
486 
487  // Skip the elements that were just refined
488  if (elem->refinement_flag() == Elem::JUST_REFINED) continue;
489 
490  for (unsigned int n=0; n<elem->n_nodes(); n++)
491  {
492  Node* node = elem->get_node(n);
493 
494  if (node->old_dof_object == NULL)
495  if (node->has_dofs(sys_num))
496  node->set_old_dof_object();
497  }
498 
499  libmesh_assert (!elem->old_dof_object);
500 
501  if (elem->has_dofs(sys_num))
502  elem->set_old_dof_object();
503  }
504  }
505 
506 #endif // #ifdef LIBMESH_ENABLE_AMR
507 
508 
509  //------------------------------------------------------------
510  // Then set the number of variables for each \p DofObject
511  // equal to n_variables() for this system. This will
512  // handle new \p DofObjects that may have just been created
513  {
514  // All the nodes
515  MeshBase::node_iterator node_it = mesh.nodes_begin();
516  const MeshBase::node_iterator node_end = mesh.nodes_end();
517 
518  for ( ; node_it != node_end; ++node_it)
519  (*node_it)->set_n_vars_per_group(sys_num, n_vars_per_group);
520 
521  // All the elements
522  MeshBase::element_iterator elem_it = mesh.elements_begin();
523  const MeshBase::element_iterator elem_end = mesh.elements_end();
524 
525  for ( ; elem_it != elem_end; ++elem_it)
526  (*elem_it)->set_n_vars_per_group(sys_num, n_vars_per_group);
527  }
528 
529 
530  // Zero _n_SCALAR_dofs, it will be updated below.
531  this->_n_SCALAR_dofs = 0;
532 
533  //------------------------------------------------------------
534  // Next allocate space for the DOF indices
535  for (unsigned int vg=0; vg<n_var_groups; vg++)
536  {
537  const VariableGroup &vg_description = this->variable_group(vg);
538 
539  const unsigned int n_var_in_group = vg_description.n_variables();
540  const FEType& base_fe_type = vg_description.type();
541 
542  // Don't need to loop over elements for a SCALAR variable
543  // Just increment _n_SCALAR_dofs
544  if(base_fe_type.family == SCALAR)
545  {
546  this->_n_SCALAR_dofs += base_fe_type.order*n_var_in_group;
547  continue;
548  }
549 
550  // This should be constant even on p-refined elements
551  const bool extra_hanging_dofs =
552  FEInterface::extra_hanging_dofs(base_fe_type);
553 
554  // For all the active elements
555  MeshBase::element_iterator elem_it = mesh.active_elements_begin();
556  const MeshBase::element_iterator elem_end = mesh.active_elements_end();
557 
558  // Count vertex degrees of freedom first
559  for ( ; elem_it != elem_end; ++elem_it)
560  {
561  Elem* elem = *elem_it;
562  libmesh_assert(elem);
563 
564  // Skip the numbering if this variable is
565  // not active on this element's subdomain
566  if (!vg_description.active_on_subdomain(elem->subdomain_id()))
567  continue;
568 
569  const ElemType type = elem->type();
570  const unsigned int dim = elem->dim();
571 
572  FEType fe_type = base_fe_type;
573 
574 #ifdef LIBMESH_ENABLE_AMR
575  // Make sure we haven't done more p refinement than we can
576  // handle
577  if (elem->p_level() + base_fe_type.order >
578  FEInterface::max_order(base_fe_type, type))
579  {
580 # ifdef DEBUG
581  if (FEInterface::max_order(base_fe_type,type) <
582  static_cast<unsigned int>(base_fe_type.order))
583  {
585  << "ERROR: Finite element "
586  << Utility::enum_to_string(base_fe_type.family)
587  << " on geometric element "
588  << Utility::enum_to_string(type) << std::endl
589  << "only supports FEInterface::max_order = "
590  << FEInterface::max_order(base_fe_type,type)
591  << ", not fe_type.order = " << base_fe_type.order
592  << std::endl;
593 
594  libmesh_error();
595  }
596 
598  << "WARNING: Finite element "
599  << Utility::enum_to_string(base_fe_type.family)
600  << " on geometric element "
601  << Utility::enum_to_string(type) << std::endl
602  << "could not be p refined past FEInterface::max_order = "
603  << FEInterface::max_order(base_fe_type,type)
604  << std::endl;
605 # endif
606  elem->set_p_level(FEInterface::max_order(base_fe_type,type)
607  - base_fe_type.order);
608  }
609 #endif
610 
611  fe_type.order = static_cast<Order>(fe_type.order +
612  elem->p_level());
613 
614  // Allocate the vertex DOFs
615  for (unsigned int n=0; n<elem->n_nodes(); n++)
616  {
617  Node* node = elem->get_node(n);
618 
619  if (elem->is_vertex(n))
620  {
621  const unsigned int old_node_dofs =
622  node->n_comp_group(sys_num, vg);
623 
624  const unsigned int vertex_dofs =
626  type, n),
627  old_node_dofs);
628 
629  // Some discontinuous FEs have no vertex dofs
630  if (vertex_dofs > old_node_dofs)
631  {
632  node->set_n_comp_group(sys_num, vg,
633  vertex_dofs);
634 
635  // Abusing dof_number to set a "this is a
636  // vertex" flag
637  node->set_vg_dof_base(sys_num, vg,
638  vertex_dofs);
639 
640  // libMesh::out << "sys_num,vg,old_node_dofs,vertex_dofs="
641  // << sys_num << ","
642  // << vg << ","
643  // << old_node_dofs << ","
644  // << vertex_dofs << '\n',
645  // node->debug_buffer();
646 
647  // libmesh_assert_equal_to (vertex_dofs, node->n_comp(sys_num, vg));
648  // libmesh_assert_equal_to (vertex_dofs, node->vg_dof_base(sys_num, vg));
649  }
650  }
651  }
652  } // done counting vertex dofs
653 
654  // count edge & face dofs next
655  elem_it = mesh.active_elements_begin();
656 
657  for ( ; elem_it != elem_end; ++elem_it)
658  {
659  Elem* elem = *elem_it;
660  libmesh_assert(elem);
661 
662  // Skip the numbering if this variable is
663  // not active on this element's subdomain
664  if (!vg_description.active_on_subdomain(elem->subdomain_id()))
665  continue;
666 
667  const ElemType type = elem->type();
668  const unsigned int dim = elem->dim();
669 
670  FEType fe_type = base_fe_type;
671  fe_type.order = static_cast<Order>(fe_type.order +
672  elem->p_level());
673 
674  // Allocate the edge and face DOFs
675  for (unsigned int n=0; n<elem->n_nodes(); n++)
676  {
677  Node* node = elem->get_node(n);
678 
679  const unsigned int old_node_dofs =
680  node->n_comp_group(sys_num, vg);
681 
682  const unsigned int vertex_dofs = old_node_dofs?
683  libmesh_cast_int<unsigned int>(node->vg_dof_base (sys_num,vg)):0;
684 
685  const unsigned int new_node_dofs =
686  FEInterface::n_dofs_at_node(dim, fe_type, type, n);
687 
688  // We've already allocated vertex DOFs
689  if (elem->is_vertex(n))
690  {
691  libmesh_assert_greater_equal (old_node_dofs, vertex_dofs);
692  // //if (vertex_dofs < new_node_dofs)
693  // libMesh::out << "sys_num,vg,old_node_dofs,vertex_dofs,new_node_dofs="
694  // << sys_num << ","
695  // << vg << ","
696  // << old_node_dofs << ","
697  // << vertex_dofs << ","
698  // << new_node_dofs << '\n',
699  // node->debug_buffer();
700 
701  libmesh_assert_greater_equal (vertex_dofs, new_node_dofs);
702  }
703  // We need to allocate the rest
704  else
705  {
706  // If this has no dofs yet, it needs no vertex
707  // dofs, so we just give it edge or face dofs
708  if (!old_node_dofs)
709  {
710  node->set_n_comp_group(sys_num, vg,
711  new_node_dofs);
712  // Abusing dof_number to set a "this has no
713  // vertex dofs" flag
714  if (new_node_dofs)
715  node->set_vg_dof_base(sys_num, vg,
716  0);
717  }
718 
719  // If this has dofs, but has no vertex dofs,
720  // it may still need more edge or face dofs if
721  // we're p-refined.
722  else if (vertex_dofs == 0)
723  {
724  if (new_node_dofs > old_node_dofs)
725  {
726  node->set_n_comp_group(sys_num, vg,
727  new_node_dofs);
728 
729  node->set_vg_dof_base(sys_num, vg,
730  vertex_dofs);
731  }
732  }
733  // If this is another element's vertex,
734  // add more (non-overlapping) edge/face dofs if
735  // necessary
736  else if (extra_hanging_dofs)
737  {
738  if (new_node_dofs > old_node_dofs - vertex_dofs)
739  {
740  node->set_n_comp_group(sys_num, vg,
741  vertex_dofs + new_node_dofs);
742 
743  node->set_vg_dof_base(sys_num, vg,
744  vertex_dofs);
745  }
746  }
747  // If this is another element's vertex, add any
748  // (overlapping) edge/face dofs if necessary
749  else
750  {
751  libmesh_assert_greater_equal (old_node_dofs, vertex_dofs);
752  if (new_node_dofs > old_node_dofs)
753  {
754  node->set_n_comp_group(sys_num, vg,
755  new_node_dofs);
756 
757  node->set_vg_dof_base (sys_num, vg,
758  vertex_dofs);
759  }
760  }
761  }
762  }
763  // Allocate the element DOFs
764  const unsigned int dofs_per_elem =
765  FEInterface::n_dofs_per_elem(dim, fe_type,
766  type);
767 
768  elem->set_n_comp_group(sys_num, vg, dofs_per_elem);
769 
770  }
771  } // end loop over variable groups
772 
773  // Calling DofMap::reinit() by itself makes little sense,
774  // so we won't bother with nonlocal DofObjects.
775  // Those will be fixed by distribute_dofs
776 
777  //------------------------------------------------------------
778  // Finally, clear all the current DOF indices
779  // (distribute_dofs expects them cleared!)
780  this->invalidate_dofs(mesh);
781 
782  STOP_LOG("reinit()", "DofMap");
783 }
void libMesh::DofMap::remove_adjoint_dirichlet_boundary ( const DirichletBoundary dirichlet_boundary,
unsigned int  q 
)

Removes from the system the specified Dirichlet boundary for the adjoint equation defined by Quantity of interest index q

Definition at line 3622 of file dof_map_constraints.C.

References libMesh::DirichletBoundary::b, end, libMesh::libmesh_assert(), libMesh::libmesh_assert_greater(), and libMesh::DirichletBoundary::variables.

3624 {
3626  qoi_index);
3627 
3628  // Find a boundary condition matching the one to be removed
3629  std::vector<DirichletBoundary *>::iterator it =
3630  _adjoint_dirichlet_boundaries[qoi_index]->begin();
3631  std::vector<DirichletBoundary *>::iterator end =
3632  _adjoint_dirichlet_boundaries[qoi_index]->end();
3633  for (; it != end; ++it)
3634  {
3635  DirichletBoundary *bdy = *it;
3636 
3637  if ((bdy->b == boundary_to_remove.b) &&
3638  bdy->variables == boundary_to_remove.variables)
3639  break;
3640  }
3641 
3642  // Delete it and remove it
3643  libmesh_assert (it != end);
3644  delete *it;
3645  _adjoint_dirichlet_boundaries[qoi_index]->erase(it);
3646 }
void libMesh::DofMap::remove_dirichlet_boundary ( const DirichletBoundary dirichlet_boundary)

Removes the specified Dirichlet boundary from the system.

Definition at line 3601 of file dof_map_constraints.C.

References libMesh::DirichletBoundary::b, end, libMesh::libmesh_assert(), and libMesh::DirichletBoundary::variables.

3602 {
3603  // Find a boundary condition matching the one to be removed
3604  std::vector<DirichletBoundary *>::iterator it = _dirichlet_boundaries->begin();
3605  std::vector<DirichletBoundary *>::iterator end = _dirichlet_boundaries->end();
3606  for (; it != end; ++it)
3607  {
3608  DirichletBoundary *bdy = *it;
3609 
3610  if ((bdy->b == boundary_to_remove.b) &&
3611  bdy->variables == boundary_to_remove.variables)
3612  break;
3613  }
3614 
3615  // Delete it and remove it
3616  libmesh_assert (it != end);
3617  delete *it;
3618  _dirichlet_boundaries->erase(it);
3619 }
void libMesh::DofMap::SCALAR_dof_indices ( std::vector< dof_id_type > &  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 1798 of file dof_map.C.

References libMesh::err, libMesh::FEType::family, libMesh::libmesh_assert(), n_dofs(), n_old_dofs(), n_SCALAR_dofs(), n_variables(), libMesh::FEType::order, libMeshEnums::SCALAR, libMesh::START_LOG(), libMesh::STOP_LOG(), libMesh::Variable::type(), and variable().

Referenced by dof_indices(), old_dof_indices(), libMesh::System::project_vector(), libMesh::System::read_parallel_data(), libMesh::System::read_SCALAR_dofs(), libMesh::System::write_parallel_data(), and libMesh::System::write_SCALAR_dofs().

1806 {
1807  START_LOG("SCALAR_dof_indices()", "DofMap");
1808 
1809  if(this->variable(vn).type().family != SCALAR)
1810  {
1811  libMesh::err << "ERROR: SCALAR_dof_indices called for a non-SCALAR variable."
1812  << std::endl;
1813  }
1814 
1815  // Clear the DOF indices vector
1816  di.clear();
1817 
1818  // First we need to find out the first dof
1819  // index for each SCALAR.
1820 #ifdef LIBMESH_ENABLE_AMR
1821  dof_id_type first_SCALAR_dof_index = (old_dofs ? n_old_dofs() : n_dofs()) - n_SCALAR_dofs();
1822 #else
1823  dof_id_type first_SCALAR_dof_index = n_dofs() - n_SCALAR_dofs();
1824 #endif
1825  std::map<unsigned int, dof_id_type> SCALAR_first_dof_index;
1826  SCALAR_first_dof_index.clear();
1827 
1828  // Iterate over _all_ of the SCALARs and store each one's first dof index
1829  // We need to do this since the SCALAR dofs are packed contiguously
1830  for (unsigned int v=0; v<this->n_variables(); v++)
1831  if(this->variable(v).type().family == SCALAR)
1832  {
1833  unsigned int current_n_SCALAR_dofs = this->variable(v).type().order;
1834  SCALAR_first_dof_index.insert(
1835  std::pair<unsigned int, dof_id_type>(v,first_SCALAR_dof_index) );
1836  first_SCALAR_dof_index += current_n_SCALAR_dofs;
1837  }
1838 
1839  // Now use vn to index into SCALAR_first_dof_index
1840  std::map<unsigned int, dof_id_type>::const_iterator iter =
1841  SCALAR_first_dof_index.find(vn);
1842