dof_map.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // Local Includes
23 #include "libmesh/dense_matrix.h"
26 #include "libmesh/dof_map.h"
27 #include "libmesh/elem.h"
28 #include "libmesh/fe_interface.h"
29 #include "libmesh/fe_type.h"
30 #include "libmesh/fe_base.h" // FEBase::build() for continuity test
33 #include "libmesh/mesh_base.h"
34 #include "libmesh/mesh_tools.h"
35 #include "libmesh/numeric_vector.h"
36 #include "libmesh/parallel.h"
38 #include "libmesh/sparse_matrix.h"
40 #include "libmesh/string_to_enum.h"
41 #include "libmesh/threads.h"
43 
44 // C++ Includes
45 #include <set>
46 #include <algorithm> // for std::fill, std::equal_range, std::max, std::lower_bound, etc.
47 #include <sstream>
48 
49 namespace libMesh
50 {
51 
52 // ------------------------------------------------------------
53 // DofMap member functions
54 UniquePtr<SparsityPattern::Build>
56 {
57  libmesh_assert (mesh.is_prepared());
58 
59  LOG_SCOPE("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.
83  (new SparsityPattern::Build (mesh,
84  *this,
85  this->_dof_coupling,
86  this->_coupling_functors,
87  implicit_neighbor_dofs,
89 
91  mesh.active_local_elements_end()), *sp);
92 
93  sp->parallel_sync();
94 
95 #ifndef NDEBUG
96  // Avoid declaring these variables unless asserts are enabled.
97  const processor_id_type proc_id = mesh.processor_id();
98  const dof_id_type n_dofs_on_proc = this->n_dofs_on_processor(proc_id);
99 #endif
100  libmesh_assert_equal_to (sp->sparsity_pattern.size(), n_dofs_on_proc);
101 
102  // Check to see if we have any extra stuff to add to the sparsity_pattern
104  {
106  {
107  libmesh_here();
108  libMesh::out << "WARNING: You have specified both an extra sparsity function and object.\n"
109  << " Are you sure this is what you meant to do??"
110  << std::endl;
111  }
112 
114  (sp->sparsity_pattern, sp->n_nz,
115  sp->n_oz, _extra_sparsity_context);
116  }
117 
120  (sp->sparsity_pattern, sp->n_nz, sp->n_oz);
121 
122  return UniquePtr<SparsityPattern::Build>(sp.release());
123 }
124 
125 
126 
127 DofMap::DofMap(const unsigned int number,
128  MeshBase & mesh) :
129  ParallelObject (mesh.comm()),
131  _variables(),
133  _sys_number(number),
134  _mesh(mesh),
135  _matrices(),
136  _first_df(),
137  _end_df(),
139  _send_list(),
151  _n_dfs(0),
152  _n_SCALAR_dofs(0)
153 #ifdef LIBMESH_ENABLE_AMR
154  , _n_old_dfs(0),
155  _first_old_df(),
156  _end_old_df(),
158 #endif
159 #ifdef LIBMESH_ENABLE_CONSTRAINTS
160  , _dof_constraints()
164 #endif
165 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
167 #endif
168 #ifdef LIBMESH_ENABLE_PERIODIC
170 #endif
171 #ifdef LIBMESH_ENABLE_DIRICHLET
174 #endif
177 {
178  _matrices.clear();
179 
180  _default_coupling->set_mesh(&_mesh);
181  _default_evaluating->set_mesh(&_mesh);
182  _default_evaluating->set_n_levels(1);
183 
184 #ifdef LIBMESH_ENABLE_PERIODIC
185  _default_coupling->set_periodic_boundaries(_periodic_boundaries.get());
186  _default_evaluating->set_periodic_boundaries(_periodic_boundaries.get());
187 #endif
188 
191 }
192 
193 
194 
195 // Destructor
197 {
198  this->clear();
199 
200  // clear() resets all but the default DofMap-based functors. We
201  // need to remove those from the mesh too before we die.
204 
205 #ifdef LIBMESH_ENABLE_DIRICHLET
206  for (std::size_t q = 0; q != _adjoint_dirichlet_boundaries.size(); ++q)
208 #endif
209 }
210 
211 
212 #ifdef LIBMESH_ENABLE_PERIODIC
213 
214 bool DofMap::is_periodic_boundary (const boundary_id_type boundaryid) const
215 {
216  if (_periodic_boundaries->count(boundaryid) != 0)
217  return true;
218 
219  return false;
220 }
221 
222 #endif
223 
224 
225 
226 // void DofMap::add_variable (const Variable & var)
227 // {
228 // libmesh_not_implemented();
229 // _variables.push_back (var);
230 // }
231 
232 
233 
235 {
236  _variable_groups.push_back(var_group);
237 
238  VariableGroup & new_var_group = _variable_groups.back();
239 
240  for (unsigned int var=0; var<new_var_group.n_variables(); var++)
241  _variables.push_back (new_var_group(var));
242 }
243 
244 
245 
247 {
248  parallel_object_only();
249 
250  // We shouldn't be trying to re-attach the same matrices repeatedly
251  libmesh_assert (std::find(_matrices.begin(), _matrices.end(),
252  &matrix) == _matrices.end());
253 
254  _matrices.push_back(&matrix);
255 
256  matrix.attach_dof_map (*this);
257 
258  // If we've already computed sparsity, then it's too late
259  // to wait for "compute_sparsity" to help with sparse matrix
260  // initialization, and we need to handle this matrix individually
261  bool computed_sparsity_already =
262  ((_n_nz && !_n_nz->empty()) ||
263  (_n_oz && !_n_oz->empty()));
264  this->comm().max(computed_sparsity_already);
265  if (computed_sparsity_already &&
267  {
268  // We'd better have already computed the full sparsity pattern
269  // if we need it here
271  libmesh_assert(_sp.get());
272 
273  matrix.update_sparsity_pattern (_sp->sparsity_pattern);
274  }
275 
276  if (matrix.need_full_sparsity_pattern())
278 }
279 
280 
281 
283 {
284  return (std::find(_matrices.begin(), _matrices.end(),
285  &matrix) != _matrices.end());
286 }
287 
288 
289 
291 {
292  return mesh.node_ptr(i);
293 }
294 
295 
296 
298 {
299  return mesh.elem_ptr(i);
300 }
301 
302 
303 
304 template <typename iterator_type>
305 void DofMap::set_nonlocal_dof_objects(iterator_type objects_begin,
306  iterator_type objects_end,
307  MeshBase & mesh,
308  dofobject_accessor objects)
309 {
310  // This function must be run on all processors at once
311  parallel_object_only();
312 
313  // First, iterate over local objects to find out how many
314  // are on each processor
315  std::vector<dof_id_type>
316  ghost_objects_from_proc(this->n_processors(), 0);
317 
318  iterator_type it = objects_begin;
319 
320  for (; it != objects_end; ++it)
321  {
322  DofObject * obj = *it;
323 
324  if (obj)
325  {
326  processor_id_type obj_procid = obj->processor_id();
327  // We'd better be completely partitioned by now
328  libmesh_assert_not_equal_to (obj_procid, DofObject::invalid_processor_id);
329  ghost_objects_from_proc[obj_procid]++;
330  }
331  }
332 
333  std::vector<dof_id_type> objects_on_proc(this->n_processors(), 0);
334  this->comm().allgather(ghost_objects_from_proc[this->processor_id()],
335  objects_on_proc);
336 
337 #ifdef DEBUG
338  for (processor_id_type p=0; p != this->n_processors(); ++p)
339  libmesh_assert_less_equal (ghost_objects_from_proc[p], objects_on_proc[p]);
340 #endif
341 
342  // Request sets to send to each processor
343  std::vector<std::vector<dof_id_type> >
344  requested_ids(this->n_processors());
345 
346  // We know how many of our objects live on each processor, so
347  // reserve() space for requests from each.
348  for (processor_id_type p=0; p != this->n_processors(); ++p)
349  if (p != this->processor_id())
350  requested_ids[p].reserve(ghost_objects_from_proc[p]);
351 
352  for (it = objects_begin; it != objects_end; ++it)
353  {
354  DofObject * obj = *it;
356  requested_ids[obj->processor_id()].push_back(obj->id());
357  }
358 #ifdef DEBUG
359  for (processor_id_type p=0; p != this->n_processors(); ++p)
360  libmesh_assert_equal_to (requested_ids[p].size(), ghost_objects_from_proc[p]);
361 #endif
362 
363  // Next set ghost object n_comps from other processors
364  for (processor_id_type p=1; p != this->n_processors(); ++p)
365  {
366  // Trade my requests with processor procup and procdown
367  processor_id_type procup =
368  cast_int<processor_id_type>((this->processor_id() + p) %
369  this->n_processors());
370  processor_id_type procdown =
371  cast_int<processor_id_type>((this->n_processors() +
372  this->processor_id() - p) %
373  this->n_processors());
374  std::vector<dof_id_type> request_to_fill;
375  this->comm().send_receive(procup, requested_ids[procup],
376  procdown, request_to_fill);
377 
378  // Fill those requests
379  const unsigned int
380  sys_num = this->sys_number(),
381  n_var_groups = this->n_variable_groups();
382 
383  std::vector<dof_id_type> ghost_data
384  (request_to_fill.size() * 2 * n_var_groups);
385 
386  for (std::size_t i=0; i != request_to_fill.size(); ++i)
387  {
388  DofObject * requested = (this->*objects)(mesh, request_to_fill[i]);
389  libmesh_assert(requested);
390  libmesh_assert_equal_to (requested->processor_id(), this->processor_id());
391  libmesh_assert_equal_to (requested->n_var_groups(sys_num), n_var_groups);
392  for (unsigned int vg=0; vg != n_var_groups; ++vg)
393  {
394  unsigned int n_comp_g =
395  requested->n_comp_group(sys_num, vg);
396  ghost_data[i*2*n_var_groups+vg] = n_comp_g;
397  dof_id_type my_first_dof = n_comp_g ?
398  requested->vg_dof_base(sys_num, vg) : 0;
399  libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
400  ghost_data[i*2*n_var_groups+n_var_groups+vg] = my_first_dof;
401  }
402  }
403 
404  // Trade back the results
405  std::vector<dof_id_type> filled_request;
406  this->comm().send_receive(procdown, ghost_data,
407  procup, filled_request);
408 
409  // And copy the id changes we've now been informed of
410  libmesh_assert_equal_to (filled_request.size(),
411  requested_ids[procup].size() * 2 * n_var_groups);
412  for (std::size_t i=0; i != requested_ids[procup].size(); ++i)
413  {
414  DofObject * requested = (this->*objects)(mesh, requested_ids[procup][i]);
415  libmesh_assert(requested);
416  libmesh_assert_equal_to (requested->processor_id(), procup);
417  for (unsigned int vg=0; vg != n_var_groups; ++vg)
418  {
419  unsigned int n_comp_g =
420  cast_int<unsigned int>(filled_request[i*2*n_var_groups+vg]);
421  requested->set_n_comp_group(sys_num, vg, n_comp_g);
422  if (n_comp_g)
423  {
424  dof_id_type my_first_dof =
425  filled_request[i*2*n_var_groups+n_var_groups+vg];
426  libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
427  requested->set_vg_dof_base
428  (sys_num, vg, my_first_dof);
429  }
430  }
431  }
432  }
433 
434 #ifdef DEBUG
435  // Double check for invalid dofs
436  for (it = objects_begin; it != objects_end; ++it)
437  {
438  DofObject * obj = *it;
439  libmesh_assert (obj);
440  unsigned int num_variables = obj->n_vars(this->sys_number());
441  for (unsigned int v=0; v != num_variables; ++v)
442  {
443  unsigned int n_comp =
444  obj->n_comp(this->sys_number(), v);
445  dof_id_type my_first_dof = n_comp ?
446  obj->dof_number(this->sys_number(), v, 0) : 0;
447  libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
448  }
449  }
450 #endif
451 }
452 
453 
454 
456 {
457  libmesh_assert (mesh.is_prepared());
458 
459  LOG_SCOPE("reinit()", "DofMap");
460 
461  // We ought to reconfigure our default coupling functor.
462  //
463  // The user might have removed it from our coupling functors set,
464  // but if so, who cares, this reconfiguration is cheap.
465  _default_coupling->set_dof_coupling(this->_dof_coupling);
466 
467  // By default we may want 0 or 1 levels of coupling
468  unsigned int standard_n_levels =
469  this->use_coupled_neighbor_dofs(mesh);
470  _default_coupling->set_n_levels
471  (std::max(_default_coupling->n_levels(), standard_n_levels));
472 
473  // But we *don't* want to restrict to a CouplingMatrix unless the
474  // user does so manually; the original libMesh behavior was to put
475  // ghost indices on the send_list regardless of variable.
476  //_default_evaluating->set_dof_coupling(this->_dof_coupling);
477 
478  const unsigned int
479  sys_num = this->sys_number(),
480  n_var_groups = this->n_variable_groups();
481 
482  // The DofObjects need to know how many variable groups we have, and
483  // how many variables there are in each group.
484  std::vector<unsigned int> n_vars_per_group; n_vars_per_group.reserve (n_var_groups);
485 
486  for (unsigned int vg=0; vg<n_var_groups; vg++)
487  n_vars_per_group.push_back (this->variable_group(vg).n_variables());
488 
489 #ifdef LIBMESH_ENABLE_AMR
490 
491  //------------------------------------------------------------
492  // Clear the old_dof_objects for all the nodes
493  // and elements so that we can overwrite them
494  {
495  MeshBase::node_iterator node_it = mesh.nodes_begin();
496  const MeshBase::node_iterator node_end = mesh.nodes_end();
497 
498  for ( ; node_it != node_end; ++node_it)
499  {
500  (*node_it)->clear_old_dof_object();
501  libmesh_assert (!(*node_it)->old_dof_object);
502  }
503 
505  const MeshBase::element_iterator elem_end = mesh.elements_end();
506 
507  for ( ; elem_it != elem_end; ++elem_it)
508  {
509  (*elem_it)->clear_old_dof_object();
510  libmesh_assert (!(*elem_it)->old_dof_object);
511  }
512  }
513 
514 
515  //------------------------------------------------------------
516  // Set the old_dof_objects for the elements that
517  // weren't just created, if these old dof objects
518  // had variables
519  {
521  const MeshBase::element_iterator elem_end = mesh.elements_end();
522 
523  for ( ; elem_it != elem_end; ++elem_it)
524  {
525  Elem * elem = *elem_it;
526 
527  // Skip the elements that were just refined
528  if (elem->refinement_flag() == Elem::JUST_REFINED) continue;
529 
530  for (unsigned int n=0; n<elem->n_nodes(); n++)
531  {
532  Node & node = elem->node_ref(n);
533 
534  if (node.old_dof_object == libmesh_nullptr)
535  if (node.has_dofs(sys_num))
536  node.set_old_dof_object();
537  }
538 
540 
541  if (elem->has_dofs(sys_num))
542  elem->set_old_dof_object();
543  }
544  }
545 
546 #endif // #ifdef LIBMESH_ENABLE_AMR
547 
548 
549  //------------------------------------------------------------
550  // Then set the number of variables for each \p DofObject
551  // equal to n_variables() for this system. This will
552  // handle new \p DofObjects that may have just been created
553  {
554  // All the nodes
555  MeshBase::node_iterator node_it = mesh.nodes_begin();
556  const MeshBase::node_iterator node_end = mesh.nodes_end();
557 
558  for ( ; node_it != node_end; ++node_it)
559  (*node_it)->set_n_vars_per_group(sys_num, n_vars_per_group);
560 
561  // All the elements
563  const MeshBase::element_iterator elem_end = mesh.elements_end();
564 
565  for ( ; elem_it != elem_end; ++elem_it)
566  (*elem_it)->set_n_vars_per_group(sys_num, n_vars_per_group);
567  }
568 
569 
570  // Zero _n_SCALAR_dofs, it will be updated below.
571  this->_n_SCALAR_dofs = 0;
572 
573  //------------------------------------------------------------
574  // Next allocate space for the DOF indices
575  for (unsigned int vg=0; vg<n_var_groups; vg++)
576  {
577  const VariableGroup & vg_description = this->variable_group(vg);
578 
579  const unsigned int n_var_in_group = vg_description.n_variables();
580  const FEType & base_fe_type = vg_description.type();
581 
582  // Don't need to loop over elements for a SCALAR variable
583  // Just increment _n_SCALAR_dofs
584  if (base_fe_type.family == SCALAR)
585  {
586  this->_n_SCALAR_dofs += base_fe_type.order.get_order()*n_var_in_group;
587  continue;
588  }
589 
590  // This should be constant even on p-refined elements
591  const bool extra_hanging_dofs =
592  FEInterface::extra_hanging_dofs(base_fe_type);
593 
594  // For all the active elements
596  const MeshBase::element_iterator elem_end = mesh.active_elements_end();
597 
598  // Count vertex degrees of freedom first
599  for ( ; elem_it != elem_end; ++elem_it)
600  {
601  Elem * elem = *elem_it;
602  libmesh_assert(elem);
603 
604  // Skip the numbering if this variable is
605  // not active on this element's subdomain
606  if (!vg_description.active_on_subdomain(elem->subdomain_id()))
607  continue;
608 
609  const ElemType type = elem->type();
610  const unsigned int dim = elem->dim();
611 
612  FEType fe_type = base_fe_type;
613 
614 #ifdef LIBMESH_ENABLE_AMR
615  // Make sure we haven't done more p refinement than we can
616  // handle
617  if (elem->p_level() + base_fe_type.order >
618  FEInterface::max_order(base_fe_type, type))
619  {
620  libmesh_assert_less_msg(static_cast<unsigned int>(base_fe_type.order.get_order()),
621  FEInterface::max_order(base_fe_type,type),
622  "ERROR: Finite element "
623  << Utility::enum_to_string(base_fe_type.family)
624  << " on geometric element "
625  << Utility::enum_to_string(type)
626  << "\nonly supports FEInterface::max_order = "
627  << FEInterface::max_order(base_fe_type,type)
628  << ", not fe_type.order = "
629  << base_fe_type.order);
630 
631 # ifdef DEBUG
632  libMesh::err << "WARNING: Finite element "
633  << Utility::enum_to_string(base_fe_type.family)
634  << " on geometric element "
635  << Utility::enum_to_string(type) << std::endl
636  << "could not be p refined past FEInterface::max_order = "
637  << FEInterface::max_order(base_fe_type,type)
638  << std::endl;
639 # endif
640  elem->set_p_level(FEInterface::max_order(base_fe_type,type)
641  - base_fe_type.order);
642  }
643 #endif
644 
645  fe_type.order = static_cast<Order>(fe_type.order +
646  elem->p_level());
647 
648  // Allocate the vertex DOFs
649  for (unsigned int n=0; n<elem->n_nodes(); n++)
650  {
651  Node & node = elem->node_ref(n);
652 
653  if (elem->is_vertex(n))
654  {
655  const unsigned int old_node_dofs =
656  node.n_comp_group(sys_num, vg);
657 
658  const unsigned int vertex_dofs =
660  type, n),
661  old_node_dofs);
662 
663  // Some discontinuous FEs have no vertex dofs
664  if (vertex_dofs > old_node_dofs)
665  {
666  node.set_n_comp_group(sys_num, vg,
667  vertex_dofs);
668 
669  // Abusing dof_number to set a "this is a
670  // vertex" flag
671  node.set_vg_dof_base(sys_num, vg,
672  vertex_dofs);
673 
674  // libMesh::out << "sys_num,vg,old_node_dofs,vertex_dofs="
675  // << sys_num << ","
676  // << vg << ","
677  // << old_node_dofs << ","
678  // << vertex_dofs << '\n',
679  // node.debug_buffer();
680 
681  // libmesh_assert_equal_to (vertex_dofs, node.n_comp(sys_num, vg));
682  // libmesh_assert_equal_to (vertex_dofs, node.vg_dof_base(sys_num, vg));
683  }
684  }
685  }
686  } // done counting vertex dofs
687 
688  // count edge & face dofs next
689  elem_it = mesh.active_elements_begin();
690 
691  for ( ; elem_it != elem_end; ++elem_it)
692  {
693  Elem * elem = *elem_it;
694  libmesh_assert(elem);
695 
696  // Skip the numbering if this variable is
697  // not active on this element's subdomain
698  if (!vg_description.active_on_subdomain(elem->subdomain_id()))
699  continue;
700 
701  const ElemType type = elem->type();
702  const unsigned int dim = elem->dim();
703 
704  FEType fe_type = base_fe_type;
705  fe_type.order = static_cast<Order>(fe_type.order +
706  elem->p_level());
707 
708  // Allocate the edge and face DOFs
709  for (unsigned int n=0; n<elem->n_nodes(); n++)
710  {
711  Node & node = elem->node_ref(n);
712 
713  const unsigned int old_node_dofs =
714  node.n_comp_group(sys_num, vg);
715 
716  const unsigned int vertex_dofs = old_node_dofs?
717  cast_int<unsigned int>(node.vg_dof_base (sys_num,vg)):0;
718 
719  const unsigned int new_node_dofs =
720  FEInterface::n_dofs_at_node(dim, fe_type, type, n);
721 
722  // We've already allocated vertex DOFs
723  if (elem->is_vertex(n))
724  {
725  libmesh_assert_greater_equal (old_node_dofs, vertex_dofs);
726  // //if (vertex_dofs < new_node_dofs)
727  // libMesh::out << "sys_num,vg,old_node_dofs,vertex_dofs,new_node_dofs="
728  // << sys_num << ","
729  // << vg << ","
730  // << old_node_dofs << ","
731  // << vertex_dofs << ","
732  // << new_node_dofs << '\n',
733  // node.debug_buffer();
734 
735  libmesh_assert_greater_equal (vertex_dofs, new_node_dofs);
736  }
737  // We need to allocate the rest
738  else
739  {
740  // If this has no dofs yet, it needs no vertex
741  // dofs, so we just give it edge or face dofs
742  if (!old_node_dofs)
743  {
744  node.set_n_comp_group(sys_num, vg,
745  new_node_dofs);
746  // Abusing dof_number to set a "this has no
747  // vertex dofs" flag
748  if (new_node_dofs)
749  node.set_vg_dof_base(sys_num, vg, 0);
750  }
751 
752  // If this has dofs, but has no vertex dofs,
753  // it may still need more edge or face dofs if
754  // we're p-refined.
755  else if (vertex_dofs == 0)
756  {
757  if (new_node_dofs > old_node_dofs)
758  {
759  node.set_n_comp_group(sys_num, vg,
760  new_node_dofs);
761 
762  node.set_vg_dof_base(sys_num, vg,
763  vertex_dofs);
764  }
765  }
766  // If this is another element's vertex,
767  // add more (non-overlapping) edge/face dofs if
768  // necessary
769  else if (extra_hanging_dofs)
770  {
771  if (new_node_dofs > old_node_dofs - vertex_dofs)
772  {
773  node.set_n_comp_group(sys_num, vg,
774  vertex_dofs + new_node_dofs);
775 
776  node.set_vg_dof_base(sys_num, vg,
777  vertex_dofs);
778  }
779  }
780  // If this is another element's vertex, add any
781  // (overlapping) edge/face dofs if necessary
782  else
783  {
784  libmesh_assert_greater_equal (old_node_dofs, vertex_dofs);
785  if (new_node_dofs > old_node_dofs)
786  {
787  node.set_n_comp_group(sys_num, vg,
788  new_node_dofs);
789 
790  node.set_vg_dof_base (sys_num, vg,
791  vertex_dofs);
792  }
793  }
794  }
795  }
796  // Allocate the element DOFs
797  const unsigned int dofs_per_elem =
798  FEInterface::n_dofs_per_elem(dim, fe_type,
799  type);
800 
801  elem->set_n_comp_group(sys_num, vg, dofs_per_elem);
802 
803  }
804  } // end loop over variable groups
805 
806  // Calling DofMap::reinit() by itself makes little sense,
807  // so we won't bother with nonlocal DofObjects.
808  // Those will be fixed by distribute_dofs
809 
810  //------------------------------------------------------------
811  // Finally, clear all the current DOF indices
812  // (distribute_dofs expects them cleared!)
813  this->invalidate_dofs(mesh);
814 }
815 
816 
817 
819 {
820  const unsigned int sys_num = this->sys_number();
821 
822  // All the nodes
823  MeshBase::node_iterator node_it = mesh.nodes_begin();
824  const MeshBase::node_iterator node_end = mesh.nodes_end();
825 
826  for ( ; node_it != node_end; ++node_it)
827  (*node_it)->invalidate_dofs(sys_num);
828 
829  // All the elements
831  const MeshBase::element_iterator elem_end = mesh.active_elements_end();
832 
833  for ( ; elem_it != elem_end; ++elem_it)
834  (*elem_it)->invalidate_dofs(sys_num);
835 }
836 
837 
838 
840 {
841  // we don't want to clear
842  // the coupling matrix!
843  // It should not change...
844  //_dof_coupling->clear();
845  //
846  // But it would be inconsistent to leave our coupling settings
847  // through a clear()...
848  _dof_coupling = NULL;
849 
850  // Reset ghosting functor statuses
851  {
852  std::set<GhostingFunctor *>::iterator gf_it = this->coupling_functors_begin();
853  const std::set<GhostingFunctor *>::iterator gf_end = this->coupling_functors_end();
854  for (; gf_it != gf_end; ++gf_it)
855  {
856  GhostingFunctor * gf = *gf_it;
857  libmesh_assert(gf);
859  }
860  this->_coupling_functors.clear();
861 
862  // Go back to default coupling
863 
864  _default_coupling->set_dof_coupling(this->_dof_coupling);
865  _default_coupling->set_n_levels(this->use_coupled_neighbor_dofs(this->_mesh));
866 
868  }
869 
870 
871  {
872  std::set<GhostingFunctor *>::iterator gf_it = this->algebraic_ghosting_functors_begin();
873  const std::set<GhostingFunctor *>::iterator gf_end = this->algebraic_ghosting_functors_end();
874  for (; gf_it != gf_end; ++gf_it)
875  {
876  GhostingFunctor * gf = *gf_it;
877  libmesh_assert(gf);
879  }
880  this->_algebraic_ghosting_functors.clear();
881 
882  // Go back to default send_list generation
883 
884  // _default_evaluating->set_dof_coupling(this->_dof_coupling);
885  _default_evaluating->set_n_levels(1);
887  }
888 
889  _variables.clear();
890  _variable_groups.clear();
891  _first_df.clear();
892  _end_df.clear();
893  _first_scalar_df.clear();
894  _send_list.clear();
895  this->clear_sparsity();
897 
898 #ifdef LIBMESH_ENABLE_AMR
899 
900  _dof_constraints.clear();
901  _stashed_dof_constraints.clear();
904  _n_old_dfs = 0;
905  _first_old_df.clear();
906  _end_old_df.clear();
907  _first_old_scalar_df.clear();
908 
909 #endif
910 
911  _matrices.clear();
912 
913  _n_dfs = 0;
914 }
915 
916 
917 
919 {
920  // This function must be run on all processors at once
921  parallel_object_only();
922 
923  // Log how long it takes to distribute the degrees of freedom
924  LOG_SCOPE("distribute_dofs()", "DofMap");
925 
926  libmesh_assert (mesh.is_prepared());
927 
928  const processor_id_type proc_id = this->processor_id();
929  const processor_id_type n_proc = this->n_processors();
930 
931  // libmesh_assert_greater (this->n_variables(), 0);
932  libmesh_assert_less (proc_id, n_proc);
933 
934  // re-init in case the mesh has changed
935  this->reinit(mesh);
936 
937  // By default distribute variables in a
938  // var-major fashion, but allow run-time
939  // specification
940  bool node_major_dofs = libMesh::on_command_line ("--node_major_dofs");
941 
942  // The DOF counter, will be incremented as we encounter
943  // new degrees of freedom
944  dof_id_type next_free_dof = 0;
945 
946  // Clear the send list before we rebuild it
947  _send_list.clear();
948 
949  // Set temporary DOF indices on this processor
950  if (node_major_dofs)
951  this->distribute_local_dofs_node_major (next_free_dof, mesh);
952  else
953  this->distribute_local_dofs_var_major (next_free_dof, mesh);
954 
955  // Get DOF counts on all processors
956  std::vector<dof_id_type> dofs_on_proc(n_proc, 0);
957  this->comm().allgather(next_free_dof, dofs_on_proc);
958 
959  // Resize and fill the _first_df and _end_df arrays
960 #ifdef LIBMESH_ENABLE_AMR
963 #endif
964 
965  _first_df.resize(n_proc);
966  _end_df.resize (n_proc);
967 
968  // Get DOF offsets
969  _first_df[0] = 0;
970  for (processor_id_type i=1; i < n_proc; ++i)
971  _first_df[i] = _end_df[i-1] = _first_df[i-1] + dofs_on_proc[i-1];
972  _end_df[n_proc-1] = _first_df[n_proc-1] + dofs_on_proc[n_proc-1];
973 
974  // Clear all the current DOF indices
975  // (distribute_dofs expects them cleared!)
976  this->invalidate_dofs(mesh);
977 
978  next_free_dof = _first_df[proc_id];
979 
980  // Set permanent DOF indices on this processor
981  if (node_major_dofs)
982  this->distribute_local_dofs_node_major (next_free_dof, mesh);
983  else
984  this->distribute_local_dofs_var_major (next_free_dof, mesh);
985 
986  libmesh_assert_equal_to (next_free_dof, _end_df[proc_id]);
987 
988  //------------------------------------------------------------
989  // At this point, all n_comp and dof_number values on local
990  // DofObjects should be correct, but a DistributedMesh might have
991  // incorrect values on non-local DofObjects. Let's request the
992  // correct values from each other processor.
993 
994  if (this->n_processors() > 1)
995  {
997  mesh.nodes_end(),
999 
1001  mesh.elements_end(),
1003  }
1004 
1005 #ifdef DEBUG
1006  {
1007  const unsigned int
1008  sys_num = this->sys_number();
1009 
1010  // Processors should all agree on DoF ids for the newly numbered
1011  // system.
1013 
1014  // DoF processor ids should match DofObject processor ids
1015  MeshBase::const_node_iterator node_it = mesh.nodes_begin();
1016  const MeshBase::const_node_iterator node_end = mesh.nodes_end();
1017  for ( ; node_it != node_end; ++node_it)
1018  {
1019  DofObject const * const dofobj = *node_it;
1020  const processor_id_type obj_proc_id = dofobj->processor_id();
1021 
1022  for (unsigned int v=0; v != dofobj->n_vars(sys_num); ++v)
1023  for (unsigned int c=0; c != dofobj->n_comp(sys_num,v); ++c)
1024  {
1025  const dof_id_type dofid = dofobj->dof_number(sys_num,v,c);
1026  libmesh_assert_greater_equal (dofid, this->first_dof(obj_proc_id));
1027  libmesh_assert_less (dofid, this->end_dof(obj_proc_id));
1028  }
1029  }
1030 
1032  const MeshBase::const_element_iterator elem_end = mesh.elements_end();
1033  for ( ; elem_it != elem_end; ++elem_it)
1034  {
1035  DofObject const * const dofobj = *elem_it;
1036  const processor_id_type obj_proc_id = dofobj->processor_id();
1037 
1038  for (unsigned int v=0; v != dofobj->n_vars(sys_num); ++v)
1039  for (unsigned int c=0; c != dofobj->n_comp(sys_num,v); ++c)
1040  {
1041  const dof_id_type dofid = dofobj->dof_number(sys_num,v,c);
1042  libmesh_assert_greater_equal (dofid, this->first_dof(obj_proc_id));
1043  libmesh_assert_less (dofid, this->end_dof(obj_proc_id));
1044  }
1045  }
1046  }
1047 #endif
1048 
1049  // Set the total number of degrees of freedom, then start finding
1050  // SCALAR degrees of freedom
1051 #ifdef LIBMESH_ENABLE_AMR
1052  _n_old_dfs = _n_dfs;
1054 #endif
1055  _n_dfs = _end_df[n_proc-1];
1056  _first_scalar_df.clear();
1058  dof_id_type current_SCALAR_dof_index = n_dofs() - n_SCALAR_dofs();
1059 
1060  // Calculate and cache the initial DoF indices for SCALAR variables.
1061  // This is an O(N_vars) calculation so we want to do it once per
1062  // renumbering rather than once per SCALAR_dof_indices() call
1063 
1064  for (unsigned int v=0; v<this->n_variables(); v++)
1065  if (this->variable(v).type().family == SCALAR)
1066  {
1067  _first_scalar_df[v] = current_SCALAR_dof_index;
1068  current_SCALAR_dof_index += this->variable(v).type().order.get_order();
1069  }
1070 
1071  // Allow our GhostingFunctor objects to reinit if necessary
1072  {
1073  std::set<GhostingFunctor *>::iterator gf_it = this->algebraic_ghosting_functors_begin();
1074  const std::set<GhostingFunctor *>::iterator gf_end = this->algebraic_ghosting_functors_end();
1075  for (; gf_it != gf_end; ++gf_it)
1076  {
1077  GhostingFunctor * gf = *gf_it;
1078  libmesh_assert(gf);
1079  gf->dofmap_reinit();
1080  }
1081  }
1082 
1083  {
1084  std::set<GhostingFunctor *>::iterator gf_it = this->coupling_functors_begin();
1085  const std::set<GhostingFunctor *>::iterator gf_end = this->coupling_functors_end();
1086  for (; gf_it != gf_end; ++gf_it)
1087  {
1088  GhostingFunctor * gf = *gf_it;
1089  libmesh_assert(gf);
1090  gf->dofmap_reinit();
1091  }
1092  }
1093 
1094  // Note that in the add_neighbors_to_send_list nodes on processor
1095  // boundaries that are shared by multiple elements are added for
1096  // each element.
1097  this->add_neighbors_to_send_list(mesh);
1098 
1099  // Here we used to clean up that data structure; now System and
1100  // EquationSystems call that for us, after we've added constraint
1101  // dependencies to the send_list too.
1102  // this->sort_send_list ();
1103 }
1104 
1105 
1106 void DofMap::local_variable_indices(std::vector<dof_id_type> & idx,
1107  const MeshBase & mesh,
1108  unsigned int var_num) const
1109 {
1110  const unsigned int sys_num = this->sys_number();
1111 
1112  // If this isn't a SCALAR variable, we need to find all its field
1113  // dofs on the mesh
1114  if (this->variable_type(var_num).family != SCALAR)
1115  {
1118 
1119  for ( ; elem_it != elem_end; ++elem_it)
1120  {
1121  // Only count dofs connected to active
1122  // elements on this processor.
1123  Elem * elem = *elem_it;
1124  const unsigned int n_nodes = elem->n_nodes();
1125 
1126  // First get any new nodal DOFS
1127  for (unsigned int n=0; n<n_nodes; n++)
1128  {
1129  Node & node = elem->node_ref(n);
1130 
1131  if (node.processor_id() < this->processor_id())
1132  continue;
1133 
1134  const unsigned int n_comp = node.n_comp(sys_num, var_num);
1135  for(unsigned int i=0; i<n_comp; i++)
1136  {
1137  const dof_id_type index = node.dof_number(sys_num,var_num,i);
1138  libmesh_assert_greater_equal (index, this->first_dof());
1139  libmesh_assert_less (index, this->end_dof());
1140 
1141  if (idx.empty() || index > idx.back())
1142  idx.push_back(index);
1143  }
1144  }
1145 
1146  // Next get any new element DOFS
1147  const unsigned int n_comp = elem->n_comp(sys_num, var_num);
1148  for (unsigned int i=0; i<n_comp; i++)
1149  {
1150  const dof_id_type index = elem->dof_number(sys_num,var_num,i);
1151  if (idx.empty() || index > idx.back())
1152  idx.push_back(index);
1153  }
1154  } // done looping over elements
1155 
1156 
1157  // we may have missed assigning DOFs to nodes that we own
1158  // but to which we have no connected elements matching our
1159  // variable restriction criterion. this will happen, for example,
1160  // if variable V is restricted to subdomain S. We may not own
1161  // any elements which live in S, but we may own nodes which are
1162  // *connected* to elements which do. in this scenario these nodes
1163  // will presently have unnumbered DOFs. we need to take care of
1164  // them here since we own them and no other processor will touch them.
1165  {
1167  const MeshBase::const_node_iterator node_end = mesh.local_nodes_end();
1168 
1169  for (; node_it != node_end; ++node_it)
1170  {
1171  Node * node = *node_it;
1172  libmesh_assert(node);
1173 
1174  const unsigned int n_comp = node->n_comp(sys_num, var_num);
1175  for (unsigned int i=0; i<n_comp; i++)
1176  {
1177  const dof_id_type index = node->dof_number(sys_num,var_num,i);
1178  if (idx.empty() || index > idx.back())
1179  idx.push_back(index);
1180  }
1181  }
1182  }
1183  }
1184  // Otherwise, count up the SCALAR dofs, if we're on the processor
1185  // that holds this SCALAR variable
1186  else if (this->processor_id() == (this->n_processors()-1))
1187  {
1188  std::vector<dof_id_type> di_scalar;
1189  this->SCALAR_dof_indices(di_scalar,var_num);
1190  idx.insert( idx.end(), di_scalar.begin(), di_scalar.end());
1191  }
1192 }
1193 
1194 
1196  MeshBase & mesh)
1197 {
1198  const unsigned int sys_num = this->sys_number();
1199  const unsigned int n_var_groups = this->n_variable_groups();
1200 
1201  //-------------------------------------------------------------------------
1202  // First count and assign temporary numbers to local dofs
1204  const MeshBase::element_iterator elem_end = mesh.active_local_elements_end();
1205 
1206  for ( ; elem_it != elem_end; ++elem_it)
1207  {
1208  // Only number dofs connected to active
1209  // elements on this processor.
1210  Elem * elem = *elem_it;
1211  const unsigned int n_nodes = elem->n_nodes();
1212 
1213  // First number the nodal DOFS
1214  for (unsigned int n=0; n<n_nodes; n++)
1215  {
1216  Node & node = elem->node_ref(n);
1217 
1218  for (unsigned vg=0; vg<n_var_groups; vg++)
1219  {
1220  const VariableGroup & vg_description(this->variable_group(vg));
1221 
1222  if ((vg_description.type().family != SCALAR) &&
1223  (vg_description.active_on_subdomain(elem->subdomain_id())))
1224  {
1225  // assign dof numbers (all at once) if this is
1226  // our node and if they aren't already there
1227  if ((node.n_comp_group(sys_num,vg) > 0) &&
1228  (node.processor_id() == this->processor_id()) &&
1229  (node.vg_dof_base(sys_num,vg) ==
1231  {
1232  node.set_vg_dof_base(sys_num, vg,
1233  next_free_dof);
1234  next_free_dof += (vg_description.n_variables()*
1235  node.n_comp_group(sys_num,vg));
1236  //node.debug_buffer();
1237  }
1238  }
1239  }
1240  }
1241 
1242  // Now number the element DOFS
1243  for (unsigned vg=0; vg<n_var_groups; vg++)
1244  {
1245  const VariableGroup & vg_description(this->variable_group(vg));
1246 
1247  if ((vg_description.type().family != SCALAR) &&
1248  (vg_description.active_on_subdomain(elem->subdomain_id())))
1249  if (elem->n_comp_group(sys_num,vg) > 0)
1250  {
1251  libmesh_assert_equal_to (elem->vg_dof_base(sys_num,vg),
1253 
1254  elem->set_vg_dof_base(sys_num,
1255  vg,
1256  next_free_dof);
1257 
1258  next_free_dof += (vg_description.n_variables()*
1259  elem->n_comp(sys_num,vg));
1260  }
1261  }
1262  } // done looping over elements
1263 
1264 
1265  // we may have missed assigning DOFs to nodes that we own
1266  // but to which we have no connected elements matching our
1267  // variable restriction criterion. this will happen, for example,
1268  // if variable V is restricted to subdomain S. We may not own
1269  // any elements which live in S, but we may own nodes which are
1270  // *connected* to elements which do. in this scenario these nodes
1271  // will presently have unnumbered DOFs. we need to take care of
1272  // them here since we own them and no other processor will touch them.
1273  {
1274  MeshBase::node_iterator node_it = mesh.local_nodes_begin();
1275  const MeshBase::node_iterator node_end = mesh.local_nodes_end();
1276 
1277  for (; node_it != node_end; ++node_it)
1278  {
1279  Node * node = *node_it;
1280  libmesh_assert(node);
1281 
1282  for (unsigned vg=0; vg<n_var_groups; vg++)
1283  {
1284  const VariableGroup & vg_description(this->variable_group(vg));
1285 
1286  if (node->n_comp_group(sys_num,vg))
1287  if (node->vg_dof_base(sys_num,vg) == DofObject::invalid_id)
1288  {
1289  node->set_vg_dof_base (sys_num,
1290  vg,
1291  next_free_dof);
1292 
1293  next_free_dof += (vg_description.n_variables()*
1294  node->n_comp(sys_num,vg));
1295  }
1296  }
1297  }
1298  }
1299 
1300  // Finally, count up the SCALAR dofs
1301  this->_n_SCALAR_dofs = 0;
1302  for (unsigned vg=0; vg<n_var_groups; vg++)
1303  {
1304  const VariableGroup & vg_description(this->variable_group(vg));
1305 
1306  if (vg_description.type().family == SCALAR)
1307  {
1308  this->_n_SCALAR_dofs += (vg_description.n_variables()*
1309  vg_description.type().order.get_order());
1310  continue;
1311  }
1312  }
1313 
1314  // Only increment next_free_dof if we're on the processor
1315  // that holds this SCALAR variable
1316  if (this->processor_id() == (this->n_processors()-1))
1317  next_free_dof += _n_SCALAR_dofs;
1318 
1319 #ifdef DEBUG
1320  {
1321  // libMesh::out << "next_free_dof=" << next_free_dof << std::endl
1322  // << "_n_SCALAR_dofs=" << _n_SCALAR_dofs << std::endl;
1323 
1324  // Make sure we didn't miss any nodes
1325  MeshTools::libmesh_assert_valid_procids<Node>(mesh);
1326 
1327  MeshBase::node_iterator node_it = mesh.local_nodes_begin();
1328  const MeshBase::node_iterator node_end = mesh.local_nodes_end();
1329  for (; node_it != node_end; ++node_it)
1330  {
1331  Node * obj = *node_it;
1332  libmesh_assert(obj);
1333  unsigned int n_var_g = obj->n_var_groups(this->sys_number());
1334  for (unsigned int vg=0; vg != n_var_g; ++vg)
1335  {
1336  unsigned int n_comp_g =
1337  obj->n_comp_group(this->sys_number(), vg);
1338  dof_id_type my_first_dof = n_comp_g ?
1339  obj->vg_dof_base(this->sys_number(), vg) : 0;
1340  libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
1341  }
1342  }
1343  }
1344 #endif // DEBUG
1345 }
1346 
1347 
1348 
1350  MeshBase & mesh)
1351 {
1352  const unsigned int sys_num = this->sys_number();
1353  const unsigned int n_var_groups = this->n_variable_groups();
1354 
1355  //-------------------------------------------------------------------------
1356  // First count and assign temporary numbers to local dofs
1357  for (unsigned vg=0; vg<n_var_groups; vg++)
1358  {
1359  const VariableGroup & vg_description(this->variable_group(vg));
1360 
1361  const unsigned int n_vars_in_group = vg_description.n_variables();
1362 
1363  // Skip the SCALAR dofs
1364  if (vg_description.type().family == SCALAR)
1365  continue;
1366 
1368  const MeshBase::element_iterator elem_end = mesh.active_local_elements_end();
1369 
1370  for ( ; elem_it != elem_end; ++elem_it)
1371  {
1372  // Only number dofs connected to active
1373  // elements on this processor.
1374  Elem * elem = *elem_it;
1375 
1376  // ... and only variables which are active on
1377  // on this element's subdomain
1378  if (!vg_description.active_on_subdomain(elem->subdomain_id()))
1379  continue;
1380 
1381  const unsigned int n_nodes = elem->n_nodes();
1382 
1383  // First number the nodal DOFS
1384  for (unsigned int n=0; n<n_nodes; n++)
1385  {
1386  Node & node = elem->node_ref(n);
1387 
1388  // assign dof numbers (all at once) if this is
1389  // our node and if they aren't already there
1390  if ((node.n_comp_group(sys_num,vg) > 0) &&
1391  (node.processor_id() == this->processor_id()) &&
1392  (node.vg_dof_base(sys_num,vg) ==
1394  {
1395  node.set_vg_dof_base(sys_num, vg, next_free_dof);
1396 
1397  next_free_dof += (n_vars_in_group*
1398  node.n_comp_group(sys_num,vg));
1399  }
1400  }
1401 
1402  // Now number the element DOFS
1403  if (elem->n_comp_group(sys_num,vg) > 0)
1404  {
1405  libmesh_assert_equal_to (elem->vg_dof_base(sys_num,vg),
1407 
1408  elem->set_vg_dof_base(sys_num,
1409  vg,
1410  next_free_dof);
1411 
1412  next_free_dof += (n_vars_in_group*
1413  elem->n_comp_group(sys_num,vg));
1414  }
1415  } // end loop on elements
1416 
1417  // we may have missed assigning DOFs to nodes that we own
1418  // but to which we have no connected elements matching our
1419  // variable restriction criterion. this will happen, for example,
1420  // if variable V is restricted to subdomain S. We may not own
1421  // any elements which live in S, but we may own nodes which are
1422  // *connected* to elements which do. in this scenario these nodes
1423  // will presently have unnumbered DOFs. we need to take care of
1424  // them here since we own them and no other processor will touch them.
1425  {
1426  MeshBase::node_iterator node_it = mesh.local_nodes_begin();
1427  const MeshBase::node_iterator node_end = mesh.local_nodes_end();
1428 
1429  for (; node_it != node_end; ++node_it)
1430  {
1431  Node * node = *node_it;
1432  libmesh_assert(node);
1433 
1434  if (node->n_comp_group(sys_num,vg))
1435  if (node->vg_dof_base(sys_num,vg) == DofObject::invalid_id)
1436  {
1437  node->set_vg_dof_base (sys_num,
1438  vg,
1439  next_free_dof);
1440 
1441  next_free_dof += (n_vars_in_group*
1442  node->n_comp_group(sys_num,vg));
1443  }
1444  }
1445  }
1446  } // end loop on variable groups
1447 
1448  // Finally, count up the SCALAR dofs
1449  this->_n_SCALAR_dofs = 0;
1450  for (unsigned vg=0; vg<n_var_groups; vg++)
1451  {
1452  const VariableGroup & vg_description(this->variable_group(vg));
1453 
1454  if (vg_description.type().family == SCALAR)
1455  {
1456  this->_n_SCALAR_dofs += (vg_description.n_variables()*
1457  vg_description.type().order.get_order());
1458  continue;
1459  }
1460  }
1461 
1462  // Only increment next_free_dof if we're on the processor
1463  // that holds this SCALAR variable
1464  if (this->processor_id() == (this->n_processors()-1))
1465  next_free_dof += _n_SCALAR_dofs;
1466 
1467 #ifdef DEBUG
1468  {
1469  // Make sure we didn't miss any nodes
1470  MeshTools::libmesh_assert_valid_procids<Node>(mesh);
1471 
1472  MeshBase::node_iterator node_it = mesh.local_nodes_begin();
1473  const MeshBase::node_iterator node_end = mesh.local_nodes_end();
1474  for (; node_it != node_end; ++node_it)
1475  {
1476  Node * obj = *node_it;
1477  libmesh_assert(obj);
1478  unsigned int n_var_g = obj->n_var_groups(this->sys_number());
1479  for (unsigned int vg=0; vg != n_var_g; ++vg)
1480  {
1481  unsigned int n_comp_g =
1482  obj->n_comp_group(this->sys_number(), vg);
1483  dof_id_type my_first_dof = n_comp_g ?
1484  obj->vg_dof_base(this->sys_number(), vg) : 0;
1485  libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
1486  }
1487  }
1488  }
1489 #endif // DEBUG
1490 }
1491 
1492 
1493 
1494 void
1495 DofMap::
1497  std::set<CouplingMatrix *> & temporary_coupling_matrices,
1498  const std::set<GhostingFunctor *>::iterator & gf_begin,
1499  const std::set<GhostingFunctor *>::iterator & gf_end,
1500  const MeshBase::const_element_iterator & elems_begin,
1501  const MeshBase::const_element_iterator & elems_end,
1503 {
1504  std::set<GhostingFunctor *>::iterator gf_it = gf_begin;
1505  for (; gf_it != gf_end; ++gf_it)
1506  {
1507  GhostingFunctor::map_type more_elements_to_ghost;
1508 
1509  GhostingFunctor * gf = *gf_it;
1510  libmesh_assert(gf);
1511  (*gf)(elems_begin, elems_end, p, more_elements_to_ghost);
1512 
1513  GhostingFunctor::map_type::iterator metg_it = more_elements_to_ghost.begin();
1514  const GhostingFunctor::map_type::iterator metg_end = more_elements_to_ghost.end();
1515  for (; metg_it != metg_end; ++metg_it)
1516  {
1517  GhostingFunctor::map_type::iterator existing_it =
1518  elements_to_ghost.find (metg_it->first);
1519  if (existing_it == elements_to_ghost.end())
1520  elements_to_ghost.insert(*metg_it);
1521  else
1522  {
1523  if (existing_it->second)
1524  {
1525  if (metg_it->second)
1526  {
1527  // If this isn't already a temporary
1528  // then we need to make one so we'll
1529  // have a non-const matrix to merge
1530  if (temporary_coupling_matrices.empty() ||
1531  temporary_coupling_matrices.find(const_cast<CouplingMatrix *>(existing_it->second)) == temporary_coupling_matrices.end())
1532  {
1533  CouplingMatrix * cm = new CouplingMatrix(*existing_it->second);
1534  temporary_coupling_matrices.insert(cm);
1535  existing_it->second = cm;
1536  }
1537  const_cast<CouplingMatrix &>(*existing_it->second) &= *metg_it->second;
1538  }
1539  else
1540  {
1541  // Any existing_it matrix merged with a full
1542  // matrix (symbolized as nullptr) gives another
1543  // full matrix (symbolizable as nullptr).
1544 
1545  // So if existing_it->second is a temporary then
1546  // we don't need it anymore; we might as well
1547  // remove it to keep the set of temporaries
1548  // small.
1549  std::set<CouplingMatrix *>::iterator temp_it =
1550  temporary_coupling_matrices.find(const_cast<CouplingMatrix *>(existing_it->second));
1551  if (temp_it != temporary_coupling_matrices.end())
1552  temporary_coupling_matrices.erase(temp_it);
1553 
1554  existing_it->second = libmesh_nullptr;
1555  }
1556  }
1557  // else we have a nullptr already, then we have a full
1558  // coupling matrix, already, and merging with anything
1559  // else won't change that, so we're done.
1560  }
1561  }
1562  }
1563 }
1564 
1565 
1566 
1568 {
1569  LOG_SCOPE("add_neighbors_to_send_list()", "DofMap");
1570 
1571  const unsigned int n_var = this->n_variables();
1572 
1573  MeshBase::const_element_iterator local_elem_it
1574  = mesh.active_local_elements_begin();
1575  const MeshBase::const_element_iterator local_elem_end
1576  = mesh.active_local_elements_end();
1577 
1578  GhostingFunctor::map_type elements_to_send;
1579 
1580  // Man, I wish we had guaranteed unique_ptr availability...
1581  std::set<CouplingMatrix *> temporary_coupling_matrices;
1582 
1583  // We need to add dofs to the send list if they've been directly
1584  // requested by an algebraic ghosting functor or they've been
1585  // indirectly requested by a coupling functor.
1586  this->merge_ghost_functor_outputs(elements_to_send,
1587  temporary_coupling_matrices,
1590  local_elem_it, local_elem_end, mesh.processor_id());
1591 
1592  this->merge_ghost_functor_outputs(elements_to_send,
1593  temporary_coupling_matrices,
1594  this->coupling_functors_begin(),
1595  this->coupling_functors_end(),
1596  local_elem_it, local_elem_end, mesh.processor_id());
1597 
1598  // Making a list of non-zero coupling matrix columns is an
1599  // O(N_var^2) operation. We cache it so we only have to do it once
1600  // per CouplingMatrix and not once per element.
1601  std::map<const CouplingMatrix *, std::vector<unsigned int> >
1602  column_variable_lists;
1603 
1604  GhostingFunctor::map_type::iterator etg_it = elements_to_send.begin();
1605  const GhostingFunctor::map_type::iterator etg_end = elements_to_send.end();
1606  for (; etg_it != etg_end; ++etg_it)
1607  {
1608  const Elem * const partner = etg_it->first;
1609 
1610  // We asked ghosting functors not to give us local elements
1611  libmesh_assert_not_equal_to
1612  (partner->processor_id(), this->processor_id());
1613 
1614  const CouplingMatrix * ghost_coupling = etg_it->second;
1615 
1616  // Loop over any present coupling matrix column variables if we
1617  // have a coupling matrix, or just add all variables to
1618  // send_list if not.
1619  if (ghost_coupling)
1620  {
1621  libmesh_assert_equal_to (ghost_coupling->size(), n_var);
1622 
1623  // Try to find a cached list of column variables.
1624  std::map<const CouplingMatrix *, std::vector<unsigned int> >::const_iterator
1625  column_variable_list = column_variable_lists.find(ghost_coupling);
1626 
1627  // If we didn't find it, then we need to create it.
1628  if (column_variable_list == column_variable_lists.end())
1629  {
1630  std::pair<std::map<const CouplingMatrix *, std::vector<unsigned int> >::iterator, bool>
1631  inserted_variable_list_pair = column_variable_lists.insert(std::make_pair(ghost_coupling,
1632  std::vector<unsigned int>()));
1633  column_variable_list = inserted_variable_list_pair.first;
1634 
1635  std::vector<unsigned int> & new_variable_list =
1636  inserted_variable_list_pair.first->second;
1637 
1638  std::vector<unsigned char> has_variable(n_var, false);
1639 
1640  for (unsigned int vi = 0; vi != n_var; ++vi)
1641  {
1642  ConstCouplingRow ccr(vi, *ghost_coupling);
1643 
1645  it = ccr.begin(),
1646  end = ccr.end();
1647  it != end; ++it)
1648  {
1649  const unsigned int vj = *it;
1650  has_variable[vj] = true;
1651  }
1652  }
1653  for (unsigned int vj = 0; vj != n_var; ++vj)
1654  {
1655  if (has_variable[vj])
1656  new_variable_list.push_back(vj);
1657  }
1658  }
1659 
1660  const std::vector<unsigned int> & variable_list =
1661  column_variable_list->second;
1662 
1663  for (std::size_t j=0; j != variable_list.size(); ++j)
1664  {
1665  const unsigned int vj=variable_list[j];
1666 
1667  std::vector<dof_id_type> di;
1668  this->dof_indices (partner, di, vj);
1669 
1670  // Insert the remote DOF indices into the send list
1671  for (std::size_t j=0; j != di.size(); ++j)
1672  if (di[j] < this->first_dof() ||
1673  di[j] >= this->end_dof())
1674  _send_list.push_back(di[j]);
1675  }
1676  }
1677  else
1678  {
1679  std::vector<dof_id_type> di;
1680  this->dof_indices (partner, di);
1681 
1682  // Insert the remote DOF indices into the send list
1683  for (std::size_t j=0; j != di.size(); ++j)
1684  if (di[j] < this->first_dof() ||
1685  di[j] >= this->end_dof())
1686  _send_list.push_back(di[j]);
1687  }
1688 
1689  }
1690 
1691  // We're now done with any merged coupling matrices we had to create.
1692  for (std::set<CouplingMatrix *>::iterator
1693  it = temporary_coupling_matrices.begin(),
1694  end = temporary_coupling_matrices.begin();
1695  it != end; ++it)
1696  delete *it;
1697 
1698  //-------------------------------------------------------------------------
1699  // Our coupling functors added dofs from neighboring elements to the
1700  // send list, but we may still need to add non-local dofs from local
1701  // elements.
1702  //-------------------------------------------------------------------------
1703 
1704  // Loop over the active local elements, adding all active elements
1705  // that neighbor an active local element to the send list.
1706  for ( ; local_elem_it != local_elem_end; ++local_elem_it)
1707  {
1708  const Elem * elem = *local_elem_it;
1709 
1710  std::vector<dof_id_type> di;
1711  this->dof_indices (elem, di);
1712 
1713  // Insert the remote DOF indices into the send list
1714  for (std::size_t j=0; j != di.size(); ++j)
1715  if (di[j] < this->first_dof() ||
1716  di[j] >= this->end_dof())
1717  _send_list.push_back(di[j]);
1718  }
1719 }
1720 
1721 
1722 
1724 {
1725  LOG_SCOPE("prepare_send_list()", "DofMap");
1726 
1727  // Check to see if we have any extra stuff to add to the send_list
1729  {
1730  if (_augment_send_list)
1731  {
1732  libmesh_here();
1733  libMesh::out << "WARNING: You have specified both an extra send list function and object.\n"
1734  << " Are you sure this is what you meant to do??"
1735  << std::endl;
1736  }
1737 
1739  }
1740 
1741  if (_augment_send_list)
1743 
1744  // First sort the send list. After this
1745  // duplicated elements will be adjacent in the
1746  // vector
1747  std::sort(_send_list.begin(), _send_list.end());
1748 
1749  // Now use std::unique to remove duplicate entries
1750  std::vector<dof_id_type>::iterator new_end =
1751  std::unique (_send_list.begin(), _send_list.end());
1752 
1753  // Remove the end of the send_list. Use the "swap trick"
1754  // from Effective STL
1755  std::vector<dof_id_type> (_send_list.begin(), new_end).swap (_send_list);
1756 }
1757 
1758 
1759 void DofMap::set_implicit_neighbor_dofs(bool implicit_neighbor_dofs)
1760 {
1762  _implicit_neighbor_dofs = implicit_neighbor_dofs;
1763 }
1764 
1765 
1767 {
1768  // If we were asked on the command line, then we need to
1769  // include sensitivities between neighbor degrees of freedom
1770  bool implicit_neighbor_dofs =
1771  libMesh::on_command_line ("--implicit_neighbor_dofs");
1772 
1773  // If the user specifies --implicit_neighbor_dofs 0, then
1774  // presumably he knows what he is doing and we won't try to
1775  // automatically turn it on even when all the variables are
1776  // discontinuous.
1777  if (implicit_neighbor_dofs)
1778  {
1779  // No flag provided defaults to 'true'
1780  int flag = 1;
1781  flag = libMesh::command_line_next ("--implicit_neighbor_dofs", flag);
1782 
1783  if (!flag)
1784  {
1785  // The user said --implicit_neighbor_dofs 0, so he knows
1786  // what he is doing and really doesn't want it.
1787  return false;
1788  }
1789  }
1790 
1791  // Possibly override the commandline option, if set_implicit_neighbor_dofs
1792  // has been called.
1794  {
1795  implicit_neighbor_dofs = _implicit_neighbor_dofs;
1796 
1797  // Again, if the user explicitly says implicit_neighbor_dofs = false,
1798  // then we return here.
1799  if (!implicit_neighbor_dofs)
1800  return false;
1801  }
1802 
1803  // Look at all the variables in this system. If every one is
1804  // discontinuous then the user must be doing DG/FVM, so be nice
1805  // and force implicit_neighbor_dofs=true.
1806  {
1807  bool all_discontinuous_dofs = true;
1808 
1809  for (unsigned int var=0; var<this->n_variables(); var++)
1810  if (FEAbstract::build (mesh.mesh_dimension(),
1811  this->variable_type(var))->get_continuity() != DISCONTINUOUS)
1812  all_discontinuous_dofs = false;
1813 
1814  if (all_discontinuous_dofs)
1815  implicit_neighbor_dofs = true;
1816  }
1817 
1818  return implicit_neighbor_dofs;
1819 }
1820 
1821 
1822 
1824 {
1825  _sp = this->build_sparsity(mesh);
1826 
1827  // It is possible that some \p SparseMatrix implementations want to
1828  // see it. Let them see it before we throw it away.
1829  std::vector<SparseMatrix<Number> *>::const_iterator
1830  pos = _matrices.begin(),
1831  end = _matrices.end();
1832 
1833  // If we need the full sparsity pattern, then we share a view of its
1834  // arrays, and we pass it in to the matrices.
1836  {
1837  _n_nz = &_sp->n_nz;
1838  _n_oz = &_sp->n_oz;
1839 
1840  for (; pos != end; ++pos)
1841  (*pos)->update_sparsity_pattern (_sp->sparsity_pattern);
1842  }
1843  // If we don't need the full sparsity pattern anymore, steal the
1844  // arrays we do need and free the rest of the memory
1845  else
1846  {
1847  if (!_n_nz)
1848  _n_nz = new std::vector<dof_id_type>();
1849  _n_nz->swap(_sp->n_nz);
1850  if (!_n_oz)
1851  _n_oz = new std::vector<dof_id_type>();
1852  _n_oz->swap(_sp->n_oz);
1853 
1854  _sp.reset();
1855  }
1856 }
1857 
1858 
1859 
1861 {
1863  {
1864  libmesh_assert(_sp.get());
1865  libmesh_assert(!_n_nz || _n_nz == &_sp->n_nz);
1866  libmesh_assert(!_n_oz || _n_oz == &_sp->n_oz);
1867  _sp.reset();
1868  }
1869  else
1870  {
1871  libmesh_assert(!_sp.get());
1872  delete _n_nz;
1873  delete _n_oz;
1874  }
1877 }
1878 
1879 
1880 
1881 void
1883 {
1884  _coupling_functors.insert(&coupling_functor);
1885  _mesh.add_ghosting_functor(coupling_functor);
1886 }
1887 
1888 
1889 
1890 void
1892 {
1893  libmesh_assert(_coupling_functors.count(&coupling_functor));
1894  _coupling_functors.erase(&coupling_functor);
1895  _mesh.remove_ghosting_functor(coupling_functor);
1896 }
1897 
1898 
1899 
1900 void
1902 {
1903  _algebraic_ghosting_functors.insert(&algebraic_ghosting_functor);
1904  _mesh.add_ghosting_functor(algebraic_ghosting_functor);
1905 }
1906 
1907 
1908 
1909 void
1911 {
1912  libmesh_assert(_algebraic_ghosting_functors.count(&algebraic_ghosting_functor));
1913  _algebraic_ghosting_functors.erase(&algebraic_ghosting_functor);
1914  _mesh.remove_ghosting_functor(algebraic_ghosting_functor);
1915 }
1916 
1917 
1918 
1920  const std::vector<dof_id_type> & dof_indices_in,
1921  DenseVectorBase<Number> & Ue) const
1922 {
1923 #ifdef LIBMESH_ENABLE_AMR
1924 
1925  // Trivial mapping
1926  libmesh_assert_equal_to (dof_indices_in.size(), Ue.size());
1927  bool has_constrained_dofs = false;
1928 
1929  for (unsigned int il=0;
1930  il != cast_int<unsigned int>(dof_indices_in.size()); il++)
1931  {
1932  const dof_id_type ig = dof_indices_in[il];
1933 
1934  if (this->is_constrained_dof (ig)) has_constrained_dofs = true;
1935 
1936  libmesh_assert_less (ig, Ug.size());
1937 
1938  Ue.el(il) = Ug(ig);
1939  }
1940 
1941  // If the element has any constrained DOFs then we need
1942  // to account for them in the mapping. This will handle
1943  // the case that the input vector is not constrained.
1944  if (has_constrained_dofs)
1945  {
1946  // Copy the input DOF indices.
1947  std::vector<dof_id_type> constrained_dof_indices(dof_indices_in);
1948 
1951 
1952  this->build_constraint_matrix_and_vector (C, H, constrained_dof_indices);
1953 
1954  libmesh_assert_equal_to (dof_indices_in.size(), C.m());
1955  libmesh_assert_equal_to (constrained_dof_indices.size(), C.n());
1956 
1957  // zero-out Ue
1958  Ue.zero();
1959 
1960  // compute Ue = C Ug, with proper mapping.
1961  const unsigned int n_original_dofs =
1962  cast_int<unsigned int>(dof_indices_in.size());
1963  for (unsigned int i=0; i != n_original_dofs; i++)
1964  {
1965  Ue.el(i) = H(i);
1966 
1967  const unsigned int n_constrained =
1968  cast_int<unsigned int>(constrained_dof_indices.size());
1969  for (unsigned int j=0; j<n_constrained; j++)
1970  {
1971  const dof_id_type jg = constrained_dof_indices[j];
1972 
1973  // If Ug is a serial or ghosted vector, then this assert is
1974  // overzealous. If Ug is a parallel vector, then this assert
1975  // is redundant.
1976  // libmesh_assert ((jg >= Ug.first_local_index()) &&
1977  // (jg < Ug.last_local_index()));
1978 
1979  Ue.el(i) += C(i,j)*Ug(jg);
1980  }
1981  }
1982  }
1983 
1984 #else
1985 
1986  // Trivial mapping
1987 
1988  const unsigned int n_original_dofs =
1989  cast_int<unsigned int>(dof_indices_in.size());
1990 
1991  libmesh_assert_equal_to (n_original_dofs, Ue.size());
1992 
1993  for (unsigned int il=0; il<n_original_dofs; il++)
1994  {
1995  const dof_id_type ig = dof_indices_in[il];
1996 
1997  libmesh_assert ((ig >= Ug.first_local_index()) && (ig < Ug.last_local_index()));
1998 
1999  Ue.el(il) = Ug(ig);
2000  }
2001 
2002 #endif
2003 }
2004 
2005 void DofMap::dof_indices (const Elem * const elem,
2006  std::vector<dof_id_type> & di) const
2007 {
2008  // We now allow elem==NULL to request just SCALAR dofs
2009  // libmesh_assert(elem);
2010 
2011  // If we are asking for current indices on an element, it ought to
2012  // be an active element (or a Side proxy, which also thinks it's
2013  // active)
2014  libmesh_assert(!elem || elem->active());
2015 
2016  LOG_SCOPE("dof_indices()", "DofMap");
2017 
2018  // Clear the DOF indices vector
2019  di.clear();
2020 
2021  const unsigned int n_vars = this->n_variables();
2022 
2023 #ifdef DEBUG
2024  // Check that sizes match in DEBUG mode
2025  std::size_t tot_size = 0;
2026 #endif
2027 
2028  if (elem && elem->type() == TRI3SUBDIVISION)
2029  {
2030  // Subdivision surface FE require the 1-ring around elem
2031  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
2032 
2033  // Ghost subdivision elements have no real dofs
2034  if (!sd_elem->is_ghost())
2035  {
2036  // Determine the nodes contributing to element elem
2037  std::vector<const Node *> elem_nodes;
2038  MeshTools::Subdivision::find_one_ring(sd_elem, elem_nodes);
2039 
2040  // Get the dof numbers
2041  for (unsigned int v=0; v<n_vars; v++)
2042  {
2043  if (this->variable(v).type().family == SCALAR &&
2044  this->variable(v).active_on_subdomain(elem->subdomain_id()))
2045  {
2046 #ifdef DEBUG
2047  tot_size += this->variable(v).type().order;
2048 #endif
2049  std::vector<dof_id_type> di_new;
2050  this->SCALAR_dof_indices(di_new,v);
2051  di.insert( di.end(), di_new.begin(), di_new.end());
2052  }
2053  else
2054  _dof_indices(elem, elem->p_level(), di, v,
2055  &elem_nodes[0], elem_nodes.size()
2056 #ifdef DEBUG
2057  , tot_size
2058 #endif
2059  );
2060  }
2061  }
2062 
2063  return;
2064  }
2065 
2066  // Get the dof numbers
2067  for (unsigned int v=0; v<n_vars; v++)
2068  {
2069  const Variable & var = this->variable(v);
2070  if (var.type().family == SCALAR &&
2071  (!elem ||
2072  var.active_on_subdomain(elem->subdomain_id())))
2073  {
2074 #ifdef DEBUG
2075  tot_size += var.type().order;
2076 #endif
2077  std::vector<dof_id_type> di_new;
2078  this->SCALAR_dof_indices(di_new,v);
2079  di.insert( di.end(), di_new.begin(), di_new.end());
2080  }
2081  else if (elem)
2082  _dof_indices(elem, elem->p_level(), di, v, elem->get_nodes(),
2083  elem->n_nodes()
2084 #ifdef DEBUG
2085  , tot_size
2086 #endif
2087  );
2088  }
2089 
2090 #ifdef DEBUG
2091  libmesh_assert_equal_to (tot_size, di.size());
2092 #endif
2093 }
2094 
2095 
2096 void DofMap::dof_indices (const Elem * const elem,
2097  std::vector<dof_id_type> & di,
2098  const unsigned int vn,
2099  int p_level) const
2100 {
2101  // We now allow elem==NULL to request just SCALAR dofs
2102  // libmesh_assert(elem);
2103 
2104  LOG_SCOPE("dof_indices()", "DofMap");
2105 
2106  // Clear the DOF indices vector
2107  di.clear();
2108 
2109  // Use the default p refinement level?
2110  if (p_level == -12345)
2111  p_level = elem ? elem->p_level() : 0;
2112 
2113 #ifdef DEBUG
2114  // Check that sizes match in DEBUG mode
2115  std::size_t tot_size = 0;
2116 #endif
2117 
2118  if (elem && elem->type() == TRI3SUBDIVISION)
2119  {
2120  // Subdivision surface FE require the 1-ring around elem
2121  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
2122 
2123  // Ghost subdivision elements have no real dofs
2124  if (!sd_elem->is_ghost())
2125  {
2126  // Determine the nodes contributing to element elem
2127  std::vector<const Node *> elem_nodes;
2128  MeshTools::Subdivision::find_one_ring(sd_elem, elem_nodes);
2129 
2130  _dof_indices(elem, p_level, di, vn, &elem_nodes[0],
2131  elem_nodes.size()
2132 #ifdef DEBUG
2133  , tot_size
2134 #endif
2135  );
2136  }
2137 
2138  return;
2139  }
2140 
2141  const Variable & var = this->variable(vn);
2142 
2143  // Get the dof numbers
2144  if (var.type().family == SCALAR &&
2145  (!elem ||
2146  var.active_on_subdomain(elem->subdomain_id())))
2147  {
2148 #ifdef DEBUG
2149  tot_size += var.type().order;
2150 #endif
2151  std::vector<dof_id_type> di_new;
2152  this->SCALAR_dof_indices(di_new,vn);
2153  di.insert( di.end(), di_new.begin(), di_new.end());
2154  }
2155  else if (elem)
2156  _dof_indices(elem, p_level, di, vn, elem->get_nodes(),
2157  elem->n_nodes()
2158 #ifdef DEBUG
2159  , tot_size
2160 #endif
2161  );
2162 
2163 #ifdef DEBUG
2164  libmesh_assert_equal_to (tot_size, di.size());
2165 #endif
2166 }
2167 
2168 
2169 void DofMap::dof_indices (const Node * const node,
2170  std::vector<dof_id_type> & di) const
2171 {
2172  // We allow node==NULL to request just SCALAR dofs
2173  // libmesh_assert(elem);
2174 
2175  LOG_SCOPE("dof_indices(Node)", "DofMap");
2176 
2177  // Clear the DOF indices vector
2178  di.clear();
2179 
2180  const unsigned int n_vars = this->n_variables();
2181  const unsigned int sys_num = this->sys_number();
2182 
2183  // Get the dof numbers
2184  for (unsigned int v=0; v<n_vars; v++)
2185  {
2186  const Variable & var = this->variable(v);
2187  if (var.type().family == SCALAR)
2188  {
2189  std::vector<dof_id_type> di_new;
2190  this->SCALAR_dof_indices(di_new,v);
2191  di.insert( di.end(), di_new.begin(), di_new.end());
2192  }
2193  else
2194  {
2195  const int n_comp = node->n_comp(sys_num,v);
2196  for (int i=0; i != n_comp; ++i)
2197  {
2198  libmesh_assert_not_equal_to
2199  (node->dof_number(sys_num,v,i),
2201  di.push_back(node->dof_number(sys_num,v,i));
2202  }
2203  }
2204  }
2205 }
2206 
2207 
2208 void DofMap::dof_indices (const Node * const node,
2209  std::vector<dof_id_type> & di,
2210  const unsigned int vn) const
2211 {
2212  if (vn == libMesh::invalid_uint)
2213  {
2214  this->dof_indices(node, di);
2215  return;
2216  }
2217 
2218  // We allow node==NULL to request just SCALAR dofs
2219  // libmesh_assert(elem);
2220 
2221  LOG_SCOPE("dof_indices(Node)", "DofMap");
2222 
2223  // Clear the DOF indices vector
2224  di.clear();
2225 
2226  const unsigned int sys_num = this->sys_number();
2227 
2228  // Get the dof numbers
2229  const Variable & var = this->variable(vn);
2230  if (var.type().family == SCALAR)
2231  {
2232  std::vector<dof_id_type> di_new;
2233  this->SCALAR_dof_indices(di_new,vn);
2234  di.insert( di.end(), di_new.begin(), di_new.end());
2235  }
2236  else
2237  {
2238  const int n_comp = node->n_comp(sys_num,vn);
2239  for (int i=0; i != n_comp; ++i)
2240  {
2241  libmesh_assert_not_equal_to
2242  (node->dof_number(sys_num,vn,i),
2244  di.push_back(node->dof_number(sys_num,vn,i));
2245  }
2246  }
2247 }
2248 
2249 
2250 void DofMap::_dof_indices (const Elem * const elem,
2251  int p_level,
2252  std::vector<dof_id_type> & di,
2253  const unsigned int v,
2254  const Node * const * nodes,
2255  unsigned int n_nodes
2256 #ifdef DEBUG
2257  ,
2258  std::size_t & tot_size
2259 #endif
2260  ) const
2261 {
2262  // This internal function is only useful on valid elements
2263  libmesh_assert(elem);
2264 
2265  const Variable & var = this->variable(v);
2266 
2267  if (var.active_on_subdomain(elem->subdomain_id()))
2268  {
2269  const ElemType type = elem->type();
2270  const unsigned int sys_num = this->sys_number();
2271  const unsigned int dim = elem->dim();
2272 
2273  // Increase the polynomial order on p refined elements
2274  FEType fe_type = var.type();
2275  fe_type.order = static_cast<Order>(fe_type.order + p_level);
2276 
2277  const bool extra_hanging_dofs =
2279 
2280 #ifdef DEBUG
2281  // The number of dofs per element is non-static for subdivision FE
2282  if (fe_type.family == SUBDIVISION)
2283  tot_size += n_nodes;
2284  else
2285  tot_size += FEInterface::n_dofs(dim,fe_type,type);
2286 #endif
2287 
2288  // Get the node-based DOF numbers
2289  for (unsigned int n=0; n<n_nodes; n++)
2290  {
2291  const Node * node = nodes[n];
2292 
2293  // There is a potential problem with h refinement. Imagine a
2294  // quad9 that has a linear FE on it. Then, on the hanging side,
2295  // it can falsely identify a DOF at the mid-edge node. This is why
2296  // we call FEInterface instead of node->n_comp() directly.
2297  const unsigned int nc = FEInterface::n_dofs_at_node (dim,
2298  fe_type,
2299  type,
2300  n);
2301 
2302  // If this is a non-vertex on a hanging node with extra
2303  // degrees of freedom, we use the non-vertex dofs (which
2304  // come in reverse order starting from the end, to
2305  // simplify p refinement)
2306  if (extra_hanging_dofs && !elem->is_vertex(n))
2307  {
2308  const int dof_offset = node->n_comp(sys_num,v) - nc;
2309 
2310  // We should never have fewer dofs than necessary on a
2311  // node unless we're getting indices on a parent element,
2312  // and we should never need the indices on such a node
2313  if (dof_offset < 0)
2314  {
2315  libmesh_assert(!elem->active());
2316  di.resize(di.size() + nc, DofObject::invalid_id);
2317  }
2318  else
2319  for (int i=node->n_comp(sys_num,v)-1; i>=dof_offset; i--)
2320  {
2321  libmesh_assert_not_equal_to (node->dof_number(sys_num,v,i),
2323  di.push_back(node->dof_number(sys_num,v,i));
2324  }
2325  }
2326  // If this is a vertex or an element without extra hanging
2327  // dofs, our dofs come in forward order coming from the
2328  // beginning
2329  else
2330  for (unsigned int i=0; i<nc; i++)
2331  {
2332  libmesh_assert_not_equal_to (node->dof_number(sys_num,v,i),
2334  di.push_back(node->dof_number(sys_num,v,i));
2335  }
2336  }
2337 
2338  // If there are any element-based DOF numbers, get them
2339  const unsigned int nc = FEInterface::n_dofs_per_elem(dim,
2340  fe_type,
2341  type);
2342  // We should never have fewer dofs than necessary on an
2343  // element unless we're getting indices on a parent element,
2344  // and we should never need those indices
2345  if (nc != 0)
2346  {
2347  if (elem->n_systems() > sys_num &&
2348  nc <= elem->n_comp(sys_num,v))
2349  {
2350  for (unsigned int i=0; i<nc; i++)
2351  {
2352  libmesh_assert_not_equal_to (elem->dof_number(sys_num,v,i),
2354 
2355  di.push_back(elem->dof_number(sys_num,v,i));
2356  }
2357  }
2358  else
2359  {
2360  libmesh_assert(!elem->active() || fe_type.family == LAGRANGE || fe_type.family == SUBDIVISION);
2361  di.resize(di.size() + nc, DofObject::invalid_id);
2362  }
2363  }
2364  }
2365 }
2366 
2367 
2368 
2369 void DofMap::SCALAR_dof_indices (std::vector<dof_id_type> & di,
2370  const unsigned int vn,
2371 #ifdef LIBMESH_ENABLE_AMR
2372  const bool old_dofs
2373 #else
2374  const bool
2375 #endif
2376  ) const
2377 {
2378  LOG_SCOPE("SCALAR_dof_indices()", "DofMap");
2379 
2380  libmesh_assert(this->variable(vn).type().family == SCALAR);
2381 
2382 #ifdef LIBMESH_ENABLE_AMR
2383  // If we're asking for old dofs then we'd better have some
2384  if (old_dofs)
2385  libmesh_assert_greater_equal(n_old_dofs(), n_SCALAR_dofs());
2386 
2387  dof_id_type my_idx = old_dofs ?
2388  this->_first_old_scalar_df[vn] : this->_first_scalar_df[vn];
2389 #else
2390  dof_id_type my_idx = this->_first_scalar_df[vn];
2391 #endif
2392 
2393  libmesh_assert_not_equal_to(my_idx, DofObject::invalid_id);
2394 
2395  // The number of SCALAR dofs comes from the variable order
2396  const int n_dofs_vn = this->variable(vn).type().order.get_order();
2397 
2398  di.resize(n_dofs_vn);
2399  for (int i = 0; i != n_dofs_vn; ++i)
2400  di[i] = my_idx++;
2401 }
2402 
2403 
2404 
2405 bool DofMap::all_semilocal_indices (const std::vector<dof_id_type> & dof_indices_in) const
2406 {
2407  // We're all semilocal unless we find a counterexample
2408  for (std::size_t i=0; i != dof_indices_in.size(); ++i)
2409  {
2410  const dof_id_type di = dof_indices_in[i];
2411  // If it's not in the local indices
2412  if (di < this->first_dof() ||
2413  di >= this->end_dof())
2414  {
2415  // and if it's not in the ghost indices, then we're not
2416  // semilocal
2417  if (!std::binary_search(_send_list.begin(), _send_list.end(), di))
2418  return false;
2419  }
2420  }
2421 
2422  return true;
2423 }
2424 
2425 
2426 template <typename DofObjectSubclass>
2427 bool DofMap::is_evaluable(const DofObjectSubclass & obj,
2428  unsigned int var_num) const
2429 {
2430  // Everything is evaluable on a local object
2431  if (obj.processor_id() == this->processor_id())
2432  return true;
2433 
2434  std::vector<dof_id_type> di;
2435 
2436  if (var_num == libMesh::invalid_uint)
2437  this->dof_indices(&obj, di);
2438  else
2439  this->dof_indices(&obj, di, var_num);
2440 
2441  // We're all semilocal unless we find a counterexample
2442  for (std::size_t i=0; i != di.size(); ++i)
2443  {
2444  const dof_id_type dof_id = di[i];
2445  // If it's not in the local indices
2446  if (dof_id < this->first_dof() ||
2447  dof_id >= this->end_dof())
2448  {
2449  // and if it's not in the ghost indices, then we're not
2450  // evaluable
2451  if (!std::binary_search(_send_list.begin(), _send_list.end(), dof_id))
2452  return false;
2453  }
2454  }
2455 
2456  return true;
2457 }
2458 
2459 
2460 
2461 #ifdef LIBMESH_ENABLE_AMR
2462 
2463 void DofMap::old_dof_indices (const Elem * const elem,
2464  std::vector<dof_id_type> & di,
2465  const unsigned int vn) const
2466 {
2467  LOG_SCOPE("old_dof_indices()", "DofMap");
2468 
2469  libmesh_assert(elem);
2470 
2471  const ElemType type = elem->type();
2472  const unsigned int sys_num = this->sys_number();
2473  const unsigned int n_vars = this->n_variables();
2474  const unsigned int dim = elem->dim();
2475 
2476  // If we have dof indices stored on the elem, and there's no chance
2477  // that we only have those indices because we were just p refined,
2478  // then we should have old dof indices too.
2479  libmesh_assert(!elem->has_dofs(sys_num) ||
2481  elem->old_dof_object);
2482 
2483  // Clear the DOF indices vector.
2484  di.clear();
2485 
2486  // Determine the nodes contributing to element elem
2487  std::vector<const Node *> elem_nodes;
2488  if (elem->type() == TRI3SUBDIVISION)
2489  {
2490  // Subdivision surface FE require the 1-ring around elem
2491  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
2492  MeshTools::Subdivision::find_one_ring(sd_elem, elem_nodes);
2493  }
2494  else
2495  {
2496  // All other FE use only the nodes of elem itself
2497  elem_nodes.resize(elem->n_nodes(), libmesh_nullptr);
2498  for (unsigned int i=0; i<elem->n_nodes(); i++)
2499  elem_nodes[i] = elem->node_ptr(i);
2500  }
2501 
2502  // Get the dof numbers
2503  for (unsigned int v=0; v<n_vars; v++)
2504  if ((v == vn) || (vn == libMesh::invalid_uint))
2505  {
2506  if (this->variable(v).type().family == SCALAR &&
2507  (!elem ||
2508  this->variable(v).active_on_subdomain(elem->subdomain_id())))
2509  {
2510  // We asked for this variable, so add it to the vector.
2511  std::vector<dof_id_type> di_new;
2512  this->SCALAR_dof_indices(di_new,v,true);
2513  di.insert( di.end(), di_new.begin(), di_new.end());
2514  }
2515  else
2516  if (this->variable(v).active_on_subdomain(elem->subdomain_id()))
2517  { // Do this for all the variables if one was not specified
2518  // or just for the specified variable
2519 
2520  // Increase the polynomial order on p refined elements,
2521  // but make sure you get the right polynomial order for
2522  // the OLD degrees of freedom
2523  int p_adjustment = 0;
2524  if (elem->p_refinement_flag() == Elem::JUST_REFINED)
2525  {
2526  libmesh_assert_greater (elem->p_level(), 0);
2527  p_adjustment = -1;
2528  }
2529  else if (elem->p_refinement_flag() == Elem::JUST_COARSENED)
2530  {
2531  p_adjustment = 1;
2532  }
2533  FEType fe_type = this->variable_type(v);
2534  fe_type.order = static_cast<Order>(fe_type.order +
2535  elem->p_level() +
2536  p_adjustment);
2537 
2538  const bool extra_hanging_dofs =
2540 
2541  // Get the node-based DOF numbers
2542  for (std::size_t n=0; n<elem_nodes.size(); n++)
2543  {
2544  const Node * node = elem_nodes[n];
2545 
2546  // There is a potential problem with h refinement. Imagine a
2547  // quad9 that has a linear FE on it. Then, on the hanging side,
2548  // it can falsely identify a DOF at the mid-edge node. This is why
2549  // we call FEInterface instead of node->n_comp() directly.
2550  const unsigned int nc = FEInterface::n_dofs_at_node (dim,
2551  fe_type,
2552  type,
2553  n);
2555 
2556  // If this is a non-vertex on a hanging node with extra
2557  // degrees of freedom, we use the non-vertex dofs (which
2558  // come in reverse order starting from the end, to
2559  // simplify p refinement)
2560  if (extra_hanging_dofs && !elem->is_vertex(n))
2561  {
2562  const int dof_offset =
2563  node->old_dof_object->n_comp(sys_num,v) - nc;
2564 
2565  // We should never have fewer dofs than necessary on a
2566  // node unless we're getting indices on a parent element
2567  // or a just-coarsened element
2568  if (dof_offset < 0)
2569  {
2570  libmesh_assert(!elem->active() || elem->refinement_flag() ==
2572  di.resize(di.size() + nc, DofObject::invalid_id);
2573  }
2574  else
2575  for (int i=node->old_dof_object->n_comp(sys_num,v)-1;
2576  i>=dof_offset; i--)
2577  {
2578  libmesh_assert_not_equal_to (node->old_dof_object->dof_number(sys_num,v,i),
2580  di.push_back(node->old_dof_object->dof_number(sys_num,v,i));
2581  }
2582  }
2583  // If this is a vertex or an element without extra hanging
2584  // dofs, our dofs come in forward order coming from the
2585  // beginning
2586  else
2587  for (unsigned int i=0; i<nc; i++)
2588  {
2589  libmesh_assert_not_equal_to (node->old_dof_object->dof_number(sys_num,v,i),
2591  di.push_back(node->old_dof_object->dof_number(sys_num,v,i));
2592  }
2593  }
2594 
2595  // If there are any element-based DOF numbers, get them
2596  const unsigned int nc = FEInterface::n_dofs_per_elem(dim,
2597  fe_type,
2598  type);
2599 
2600  // We should never have fewer dofs than necessary on an
2601  // element unless we're getting indices on a parent element
2602  // or a just-coarsened element
2603  if (nc != 0)
2604  {
2605  if (elem->old_dof_object->n_systems() > sys_num &&
2606  nc <= elem->old_dof_object->n_comp(sys_num,v))
2607  {
2609 
2610  for (unsigned int i=0; i<nc; i++)
2611  {
2612  libmesh_assert_not_equal_to (elem->old_dof_object->dof_number(sys_num,v,i),
2614 
2615  di.push_back(elem->old_dof_object->dof_number(sys_num,v,i));
2616  }
2617  }
2618  else
2619  {
2620  libmesh_assert(!elem->active() || fe_type.family == LAGRANGE ||
2622  di.resize(di.size() + nc, DofObject::invalid_id);
2623  }
2624  }
2625  }
2626  } // end loop over variables
2627 }
2628 
2629 
2630 
2631 /*
2632  void DofMap::augment_send_list_for_projection(const MeshBase & mesh)
2633  {
2634  // Loop over the active local elements in the mesh.
2635  // If the element was just refined and its parent lives
2636  // on a different processor then we need to augment the
2637  // _send_list with the parent's DOF indices so that we
2638  // can access the parent's data for computing solution
2639  // projections, etc...
2640 
2641  // The DOF indices for the parent
2642  std::vector<dof_id_type> di;
2643 
2644  // Flag telling us if we need to re-sort the send_list.
2645  // Note we won't need to re-sort unless a child with a
2646  // parent belonging to a different processor is found
2647  bool needs_sorting = false;
2648 
2649 
2650  MeshBase::const_element_iterator elem_it = mesh.active_local_elements_begin();
2651  const MeshBase::const_element_iterator elem_end = mesh.active_local_elements_end();
2652 
2653  for ( ; elem_it != elem_end; ++elem_it)
2654  {
2655  const Elem * const elem = *elem_it;
2656 
2657  // We only need to consider the children that
2658  // were just refined
2659  if (elem->refinement_flag() != Elem::JUST_REFINED) continue;
2660 
2661  const Elem * const parent = elem->parent();
2662 
2663  // If the parent lives on another processor
2664  // than the child
2665  if (parent != libmesh_nullptr)
2666  if (parent->processor_id() != elem->processor_id())
2667  {
2668  // Get the DOF indices for the parent
2669  this->dof_indices (parent, di);
2670 
2671  // Insert the DOF indices into the send list
2672  _send_list.insert (_send_list.end(),
2673  di.begin(), di.end());
2674 
2675  // We will need to re-sort the send list
2676  needs_sorting = true;
2677  }
2678  }
2679 
2680  // The send-list might need to be sorted again.
2681  if (needs_sorting) this->sort_send_list ();
2682  }
2683 */
2684 
2685 #endif // LIBMESH_ENABLE_AMR
2686 
2687 
2688 #ifdef LIBMESH_ENABLE_CONSTRAINTS
2689 
2690 void DofMap::find_connected_dofs (std::vector<dof_id_type> & elem_dofs) const
2691 {
2692  typedef std::set<dof_id_type> RCSet;
2693 
2694  // First insert the DOFS we already depend on into the set.
2695  RCSet dof_set (elem_dofs.begin(), elem_dofs.end());
2696 
2697  bool done = true;
2698 
2699  // Next insert any dofs those might be constrained in terms
2700  // of. Note that in this case we may not be done: Those may
2701  // in turn depend on others. So, we need to repeat this process
2702  // in that case until the system depends only on unconstrained
2703  // degrees of freedom.
2704  for (std::size_t i=0; i<elem_dofs.size(); i++)
2705  if (this->is_constrained_dof(elem_dofs[i]))
2706  {
2707  // If the DOF is constrained
2708  DofConstraints::const_iterator
2709  pos = _dof_constraints.find(elem_dofs[i]);
2710 
2711  libmesh_assert (pos != _dof_constraints.end());
2712 
2713  const DofConstraintRow & constraint_row = pos->second;
2714 
2715  // adaptive p refinement currently gives us lots of empty constraint
2716  // rows - we should optimize those DoFs away in the future. [RHS]
2717  //libmesh_assert (!constraint_row.empty());
2718 
2719  DofConstraintRow::const_iterator it = constraint_row.begin();
2720  DofConstraintRow::const_iterator it_end = constraint_row.end();
2721 
2722 
2723  // Add the DOFs this dof is constrained in terms of.
2724  // note that these dofs might also be constrained, so
2725  // we will need to call this function recursively.
2726  for ( ; it != it_end; ++it)
2727  if (!dof_set.count (it->first))
2728  {
2729  dof_set.insert (it->first);
2730  done = false;
2731  }
2732  }
2733 
2734 
2735  // If not done then we need to do more work
2736  // (obviously :-) )!
2737  if (!done)
2738  {
2739  // Fill the vector with the contents of the set
2740  elem_dofs.clear();
2741  elem_dofs.insert (elem_dofs.end(),
2742  dof_set.begin(), dof_set.end());
2743 
2744 
2745  // May need to do this recursively. It is possible
2746  // that we just replaced a constrained DOF with another
2747  // constrained DOF.
2748  this->find_connected_dofs (elem_dofs);
2749 
2750  } // end if (!done)
2751 }
2752 
2753 #endif // LIBMESH_ENABLE_CONSTRAINTS
2754 
2755 
2756 
2757 void DofMap::print_info(std::ostream & os) const
2758 {
2759  os << this->get_info();
2760 }
2761 
2762 
2763 
2764 std::string DofMap::get_info() const
2765 {
2766  std::ostringstream os;
2767 
2768  // If we didn't calculate the exact sparsity pattern, the threaded
2769  // sparsity pattern assembly may have just given us an upper bound
2770  // on sparsity.
2771  const char * may_equal = " <= ";
2772 
2773  // If we calculated the exact sparsity pattern, then we can report
2774  // exact bandwidth figures:
2775  std::vector<SparseMatrix<Number> * >::const_iterator
2776  pos = _matrices.begin(),
2777  end = _matrices.end();
2778 
2779  for (; pos != end; ++pos)
2780  if ((*pos)->need_full_sparsity_pattern())
2781  may_equal = " = ";
2782 
2783  dof_id_type max_n_nz = 0, max_n_oz = 0;
2784  long double avg_n_nz = 0, avg_n_oz = 0;
2785 
2786  if (_n_nz)
2787  {
2788  for (std::size_t i = 0; i != _n_nz->size(); ++i)
2789  {
2790  max_n_nz = std::max(max_n_nz, (*_n_nz)[i]);
2791  avg_n_nz += (*_n_nz)[i];
2792  }
2793 
2794  std::size_t n_nz_size = _n_nz->size();
2795 
2796  this->comm().max(max_n_nz);
2797  this->comm().sum(avg_n_nz);
2798  this->comm().sum(n_nz_size);
2799 
2800  avg_n_nz /= std::max(n_nz_size,std::size_t(1));
2801 
2803 
2804  for (std::size_t i = 0; i != (*_n_oz).size(); ++i)
2805  {
2806  max_n_oz = std::max(max_n_oz, (*_n_oz)[i]);
2807  avg_n_oz += (*_n_oz)[i];
2808  }
2809 
2810  std::size_t n_oz_size = _n_oz->size();
2811 
2812  this->comm().max(max_n_oz);
2813  this->comm().sum(avg_n_oz);
2814  this->comm().sum(n_oz_size);
2815 
2816  avg_n_oz /= std::max(n_oz_size,std::size_t(1));
2817  }
2818 
2819  os << " DofMap Sparsity\n Average On-Processor Bandwidth"
2820  << may_equal << avg_n_nz << '\n';
2821 
2822  os << " Average Off-Processor Bandwidth"
2823  << may_equal << avg_n_oz << '\n';
2824 
2825  os << " Maximum On-Processor Bandwidth"
2826  << may_equal << max_n_nz << '\n';
2827 
2828  os << " Maximum Off-Processor Bandwidth"
2829  << may_equal << max_n_oz << std::endl;
2830 
2831 #ifdef LIBMESH_ENABLE_CONSTRAINTS
2832 
2833  std::size_t n_constraints = 0, max_constraint_length = 0,
2834  n_rhss = 0;
2835  long double avg_constraint_length = 0.;
2836 
2837  for (DofConstraints::const_iterator it=_dof_constraints.begin();
2838  it != _dof_constraints.end(); ++it)
2839  {
2840  // Only count local constraints, then sum later
2841  const dof_id_type constrained_dof = it->first;
2842  if (constrained_dof < this->first_dof() ||
2843  constrained_dof >= this->end_dof())
2844  continue;
2845 
2846  const DofConstraintRow & row = it->second;
2847  std::size_t rowsize = row.size();
2848 
2849  max_constraint_length = std::max(max_constraint_length,
2850  rowsize);
2851  avg_constraint_length += rowsize;
2852  n_constraints++;
2853 
2854  if (_primal_constraint_values.count(constrained_dof))
2855  n_rhss++;
2856  }
2857 
2858  this->comm().sum(n_constraints);
2859  this->comm().sum(n_rhss);
2860  this->comm().sum(avg_constraint_length);
2861  this->comm().max(max_constraint_length);
2862 
2863  os << " DofMap Constraints\n Number of DoF Constraints = "
2864  << n_constraints;
2865  if (n_rhss)
2866  os << '\n'
2867  << " Number of Heterogenous Constraints= " << n_rhss;
2868  if (n_constraints)
2869  {
2870  avg_constraint_length /= n_constraints;
2871 
2872  os << '\n'
2873  << " Average DoF Constraint Length= " << avg_constraint_length;
2874  }
2875 
2876 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2877  std::size_t n_node_constraints = 0, max_node_constraint_length = 0,
2878  n_node_rhss = 0;
2879  long double avg_node_constraint_length = 0.;
2880 
2881  for (NodeConstraints::const_iterator it=_node_constraints.begin();
2882  it != _node_constraints.end(); ++it)
2883  {
2884  // Only count local constraints, then sum later
2885  const Node * node = it->first;
2886  if (node->processor_id() != this->processor_id())
2887  continue;
2888 
2889  const NodeConstraintRow & row = it->second.first;
2890  std::size_t rowsize = row.size();
2891 
2892  max_node_constraint_length = std::max(max_node_constraint_length,
2893  rowsize);
2894  avg_node_constraint_length += rowsize;
2895  n_node_constraints++;
2896 
2897  if (it->second.second != Point(0))
2898  n_node_rhss++;
2899  }
2900 
2901  this->comm().sum(n_node_constraints);
2902  this->comm().sum(n_node_rhss);
2903  this->comm().sum(avg_node_constraint_length);
2904  this->comm().max(max_node_constraint_length);
2905 
2906  os << "\n Number of Node Constraints = " << n_node_constraints;
2907  if (n_node_rhss)
2908  os << '\n'
2909  << " Number of Heterogenous Node Constraints= " << n_node_rhss;
2910  if (n_node_constraints)
2911  {
2912  avg_node_constraint_length /= n_node_constraints;
2913  os << "\n Maximum Node Constraint Length= " << max_node_constraint_length
2914  << '\n'
2915  << " Average Node Constraint Length= " << avg_node_constraint_length;
2916  }
2917 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2918 
2919  os << std::endl;
2920 
2921 #endif // LIBMESH_ENABLE_CONSTRAINTS
2922 
2923  return os.str();
2924 }
2925 
2926 
2927 template bool DofMap::is_evaluable<Elem>(const Elem &, unsigned int) const;
2928 template bool DofMap::is_evaluable<Node>(const Node &, unsigned int) const;
2929 
2930 } // namespace libMesh
std::vector< VariableGroup > _variable_groups
Definition: dof_map.h:1392
unsigned int n_comp_group(const unsigned int s, const unsigned int vg) const
Definition: dof_object.h:785
static unsigned int n_dofs_per_elem(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:460
bool has_dofs(const unsigned int s=libMesh::invalid_uint) const
Definition: dof_object.h:847
Manages the family, order, etc. parameters for a given FE.
Definition: fe_type.h:178
DofObject * elem_ptr(MeshBase &mesh, dof_id_type i) const
Definition: dof_map.C:297
void set_p_level(const unsigned int p)
Definition: elem.h:2343
FEFamily family
Definition: fe_type.h:203
int get_order() const
Definition: fe_type.h:77
bool _implicit_neighbor_dofs_initialized
Definition: dof_map.h:1620
void distribute_local_dofs_node_major(dof_id_type &next_free_dof, MeshBase &mesh)
Definition: dof_map.C:1195
const FEType & type() const
Definition: variable.h:119
const unsigned int _sys_number
Definition: dof_map.h:1397
bool _implicit_neighbor_dofs
Definition: dof_map.h:1621
std::string get_info() const
Definition: dof_map.C:2764
void print_info(std::ostream &os=libMesh::out) const
Definition: dof_map.C:2757
UniquePtr< DefaultCoupling > _default_coupling
Definition: dof_map.h:1471
A geometric point in (x,y,z) space associated with a DOF.
Definition: node.h:52
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:414
void set_old_dof_object()
Definition: dof_object.C:142
std::set< GhostingFunctor * >::const_iterator algebraic_ghosting_functors_begin() const
Definition: dof_map.h:311
bool active() const
Definition: elem.h:2053
unsigned int n() const
virtual numeric_index_type size() const =0
const unsigned int invalid_uint
Definition: libmesh.h:184
virtual numeric_index_type last_local_index() const =0
virtual ElemType type() const =0
LIBMESH_BEST_UNORDERED_MAP< const Elem *, const CouplingMatrix * > map_type
unsigned int n_vars(const unsigned int s, const unsigned int vg) const
Definition: dof_object.h:737
void * _extra_sparsity_context
Definition: dof_map.h:1448
unsigned int p_level() const
Definition: elem.h:2218
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1676
bool is_evaluable(const DofObjectSubclass &obj, unsigned int var_num=libMesh::invalid_uint) const
Definition: dof_map.C:2427
void set_implicit_neighbor_dofs(bool implicit_neighbor_dofs)
Definition: dof_map.C:1759
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
Maps between boundary ids and PeriodicBoundaryBase objects.
std::set< GhostingFunctor * >::const_iterator coupling_functors_begin() const
Definition: dof_map.h:275
std::vector< dof_id_type > _send_list
Definition: dof_map.h:1431
bool is_attached(SparseMatrix< Number > &matrix)
Definition: dof_map.C:282
void attach_matrix(SparseMatrix< Number > &matrix)
Definition: dof_map.C:246
processor_id_type n_processors() const
void remove_ghosting_functor(GhostingFunctor &ghosting_functor)
Definition: mesh_base.C:302
DofObject * node_ptr(MeshBase &mesh, dof_id_type i) const
Definition: dof_map.C:290
The base class for all geometric element types.
Definition: elem.h:86
MeshBase & mesh
uint8_t processor_id_type
Definition: id_types.h:99
void _dof_indices(const Elem *const elem, int p_level, std::vector< dof_id_type > &di, const unsigned int v, const Node *const *nodes, unsigned int n_nodes#ifdef DEBUG, std::size_t &tot_size#endif) const
Definition: dof_map.C:2250
const class libmesh_nullptr_t libmesh_nullptr
AugmentSparsityPattern * _augment_sparsity_pattern
Definition: dof_map.h:1436
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:802
virtual const Node * node_ptr(const dof_id_type i) const =0
std::set< GhostingFunctor * > _algebraic_ghosting_functors
Definition: dof_map.h:1489
std::vector< dof_id_type > _end_old_df
Definition: dof_map.h:1556
OrderWrapper order
Definition: fe_type.h:197
void add_coupling_functor(GhostingFunctor &coupling_functor)
Definition: dof_map.C:1882
IterBase * end
void local_variable_indices(std::vector< dof_id_type > &idx, const MeshBase &mesh, unsigned int var_num) const
Definition: dof_map.C:1106
Utility class for defining generic ranges for threading.
Definition: stored_range.h:52
void set_vg_dof_base(const unsigned int s, const unsigned int vg, const dof_id_type db)
Definition: dof_object.h:894
void invalidate_dofs(MeshBase &mesh) const
Definition: dof_map.C:818
virtual void augment_send_list(std::vector< dof_id_type > &send_list)=0
void old_dof_indices(const Elem *const elem, std::vector< dof_id_type > &di, const unsigned int vn=libMesh::invalid_uint) const
Definition: dof_map.C:2463
const_iterator end() const
const Variable & variable(const unsigned int c) const
Definition: dof_map.h:1646
const unsigned int n_vars
Definition: tecplot_io.C:68
dof_id_type dof_id
Definition: xdr_io.C:48
UniquePtr< DefaultCoupling > _default_evaluating
Definition: dof_map.h:1479
dof_id_type n_dofs() const
Definition: dof_map.h:510
virtual void zero()=0
const VariableGroup & variable_group(const unsigned int c) const
Definition: dof_map.h:1636
long double max(long double a, double b)
std::vector< dof_id_type > _first_old_df
Definition: dof_map.h:1551
static unsigned int max_order(const FEType &fe_t, const ElemType &el_t)
Base class for Mesh.
Definition: mesh_base.h:67
std::vector< dof_id_type > _first_scalar_df
Definition: dof_map.h:1425
dof_id_type n_dofs_on_processor(const processor_id_type proc) const
Definition: dof_map.h:526
AdjointDofConstraintValues _adjoint_constraint_values
Definition: dof_map.h:1575
libmesh_assert(j)
std::vector< dof_id_type > * _n_nz
Definition: dof_map.h:1522
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
AugmentSendList * _augment_send_list
Definition: dof_map.h:1453
std::vector< dof_id_type > * _n_oz
Definition: dof_map.h:1528
bool need_full_sparsity_pattern
Definition: dof_map.h:1508
virtual unsigned int n_nodes() const =0
virtual node_iterator nodes_begin()=0
void add_algebraic_ghosting_functor(GhostingFunctor &ghosting_functor)
Definition: dof_map.C:1901
virtual element_iterator elements_begin()=0
dof_id_type vg_dof_base(const unsigned int s, const unsigned int vg) const
Definition: dof_object.h:914
UniquePtr< SparsityPattern::Build > _sp
Definition: dof_map.h:1514
void send_receive(const unsigned int dest_processor_id, const T1 &send, const unsigned int source_processor_id, T2 &recv, const MessageTag &send_tag=no_tag, const MessageTag &recv_tag=any_tag) const
const dof_id_type n_nodes
Definition: tecplot_io.C:67
bool is_constrained_dof(const dof_id_type dof) const
Definition: dof_map.h:1713
A variable which is solved for in a System of equations.
Definition: variable.h:49
void add_neighbors_to_send_list(MeshBase &mesh)
Definition: dof_map.C:1567
std::vector< dof_id_type > _first_old_scalar_df
Definition: dof_map.h:1562
int8_t boundary_id_type
Definition: id_types.h:51
void libmesh_assert_valid_dof_ids(const MeshBase &mesh, unsigned int sysnum)
Definition: mesh_tools.C:1379
void(* _extra_sparsity_function)(SparsityPattern::Graph &, std::vector< dof_id_type > &n_nz, std::vector< dof_id_type > &n_oz, void *)
Definition: dof_map.h:1441
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1724
dof_id_type _n_SCALAR_dofs
Definition: dof_map.h:1539
void set_nonlocal_dof_objects(iterator_type objects_begin, iterator_type objects_end, MeshBase &mesh, dofobject_accessor objects)
Definition: dof_map.C:305
virtual node_iterator local_nodes_end()=0
void remove_algebraic_ghosting_functor(GhostingFunctor &ghosting_functor)
Definition: dof_map.C:1910
static bool extra_hanging_dofs(const FEType &fe_t)
unsigned int m() const
static const processor_id_type invalid_processor_id
Definition: dof_object.h:345
virtual element_iterator active_local_elements_begin()=0
bool is_periodic_boundary(const boundary_id_type boundaryid) const
Definition: dof_map.C:214
virtual T el(const unsigned int i) const =0
unsigned int n_var_groups(const unsigned int s) const
Definition: dof_object.h:727
unsigned int n_variables() const
Definition: variable.h:217
virtual element_iterator elements_end()=0
UniquePtr< SparsityPattern::Build > build_sparsity(const MeshBase &mesh) const
Definition: dof_map.C:55
bool active_on_subdomain(subdomain_id_type sid) const
Definition: variable.h:136
virtual bool need_full_sparsity_pattern() const
void clear()
Definition: dof_map.C:839
const Node & node_ref(const unsigned int i) const
Definition: elem.h:1746
A surface shell element used in mechanics calculations.
static const dof_id_type invalid_id
Definition: dof_object.h:334
void distribute_dofs(MeshBase &)
Definition: dof_map.C:918
RefinementState p_refinement_flag() const
Definition: elem.h:2305
DofConstraints _dof_constraints
Definition: dof_map.h:1571
void reinit(MeshBase &mesh)
Definition: dof_map.C:455
DofMap(const unsigned int sys_number, MeshBase &mesh)
Definition: dof_map.C:127
CouplingMatrix * _dof_coupling
Definition: dof_map.h:1219
unsigned int n_variable_groups() const
Definition: dof_map.h:469
OStreamProxy err(std::cerr)
subdomain_id_type subdomain_id() const
Definition: elem.h:1799
std::vector< DirichletBoundaries * > _adjoint_dirichlet_boundaries
Definition: dof_map.h:1611
unsigned int sys_number() const
Definition: dof_map.h:1628
bool is_prepared() const
Definition: mesh_base.h:130
void * _extra_send_list_context
Definition: dof_map.h:1463
static void merge_ghost_functor_outputs(GhostingFunctor::map_type &elements_to_ghost, std::set< CouplingMatrix * > &temporary_coupling_matrices, const std::set< GhostingFunctor * >::iterator &gf_begin, const std::set< GhostingFunctor * >::iterator &gf_end, const MeshBase::const_element_iterator &elems_begin, const MeshBase::const_element_iterator &elems_end, processor_id_type p)
Definition: dof_map.C:1496
const_iterator begin() const
void set_n_comp_group(const unsigned int s, const unsigned int vg, const unsigned int ncomp)
Definition: dof_object.C:349
virtual element_iterator active_elements_begin()=0
std::set< GhostingFunctor * >::const_iterator algebraic_ghosting_functors_end() const
Definition: dof_map.h:317
std::string enum_to_string(const T e)
std::vector< SparseMatrix< Number > * > _matrices
Definition: dof_map.h:1409
dof_id_type end_dof() const
Definition: dof_map.h:578
virtual element_iterator active_elements_end()=0
void attach_dof_map(const DofMap &dof_map)
static unsigned int n_dofs_at_node(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int n)
Definition: fe_interface.C:436
unsigned int n_variables() const
Definition: dof_map.h:477
virtual void augment_sparsity_pattern(SparsityPattern::Graph &sparsity, std::vector< dof_id_type > &n_nz, std::vector< dof_id_type > &n_oz)=0
virtual node_iterator local_nodes_begin()=0
T command_line_next(const std::string &name, T value)
Definition: libmesh.C:957
UniquePtr< DirichletBoundaries > _dirichlet_boundaries
Definition: dof_map.h:1605
void SCALAR_dof_indices(std::vector< dof_id_type > &di, const unsigned int vn, const bool old_dofs=false) const
Definition: dof_map.C:2369
void remove_coupling_functor(GhostingFunctor &coupling_functor)
Definition: dof_map.C:1891
virtual void update_sparsity_pattern(const SparsityPattern::Graph &)
virtual node_iterator nodes_end()=0
DofObject * old_dof_object
Definition: dof_object.h:89
void find_one_ring(const Tri3Subdivision *elem, std::vector< const Node * > &nodes)
virtual bool is_vertex(const unsigned int i) const =0
dof_id_type n_old_dofs() const
Definition: dof_map.h:1172
virtual unsigned int size() const =0
X_input swap(X_system)
bool on_command_line(const std::string &arg)
Definition: libmesh.C:914
const Parallel::Communicator & comm() const
DofConstraintValueMap _primal_constraint_values
Definition: dof_map.h:1573
std::vector< Variable > _variables
Definition: dof_map.h:1387
unsigned int mesh_dimension() const
Definition: mesh_base.C:147
DofConstraints _stashed_dof_constraints
Definition: dof_map.h:1571
std::map< dof_id_type, Real, std::less< dof_id_type >, Threads::scalable_allocator< std::pair< const dof_id_type, Real > > > DofConstraintRow
Definition: dof_map.h:88
virtual unsigned int dim() const =0
dof_id_type _n_old_dfs
Definition: dof_map.h:1546
void clear_sparsity()
Definition: dof_map.C:1860
bool all_semilocal_indices(const std::vector< dof_id_type > &dof_indices) const
Definition: dof_map.C:2405
RefinementState refinement_flag() const
Definition: elem.h:2289
void parallel_reduce(const Range &range, Body &body)
Definition: threads_none.h:101
std::set< GhostingFunctor * > _coupling_functors
Definition: dof_map.h:1502
virtual numeric_index_type first_local_index() const =0
UniquePtr< PeriodicBoundaries > _periodic_boundaries
Definition: dof_map.h:1591
std::vector< dof_id_type > _end_df
Definition: dof_map.h:1419
dof_id_type n_SCALAR_dofs() const
Definition: dof_map.h:515
std::vector< dof_id_type > _first_df
Definition: dof_map.h:1414
void(* _extra_send_list_function)(std::vector< dof_id_type > &, void *)
Definition: dof_map.h:1458
DofObject *(DofMap::* dofobject_accessor)(MeshBase &mesh, dof_id_type i) const
Definition: dof_map.h:1271
std::map< const Node *, Real, std::less< const Node * >, Threads::scalable_allocator< std::pair< const Node *const, Real > > > NodeConstraintRow
Definition: dof_map.h:136
MeshBase & _mesh
Definition: dof_map.h:1402
void find_connected_dofs(std::vector< dof_id_type > &elem_dofs) const
Definition: dof_map.C:2690
void prepare_send_list()
Definition: dof_map.C:1723
Real size() const
Definition: type_vector.h:893
static UniquePtr< FEAbstract > build(const unsigned int dim, const FEType &type)
Definition: fe_abstract.C:44
std::set< GhostingFunctor * >::const_iterator coupling_functors_end() const
Definition: dof_map.h:281
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:772
void compute_sparsity(const MeshBase &)
Definition: dof_map.C:1823
dof_id_type id() const
Definition: dof_object.h:624
unsigned int n_systems() const
Definition: dof_object.h:718
virtual const Elem * elem_ptr(const dof_id_type i) const =0
OStreamProxy out(std::cout)
bool use_coupled_neighbor_dofs(const MeshBase &mesh) const
Definition: dof_map.C:1766
virtual element_iterator active_local_elements_end()=0
void extract_local_vector(const NumericVector< Number > &Ug, const std::vector< dof_id_type > &dof_indices, DenseVectorBase< Number > &Ue) const
Definition: dof_map.C:1919
void add_variable_group(const VariableGroup &var_group)
Definition: dof_map.C:234
A geometric point in (x,y,z) space.
Definition: point.h:38
const Node *const * get_nodes() const
Definition: elem.h:1716
void distribute_local_dofs_var_major(dof_id_type &next_free_dof, MeshBase &mesh)
Definition: dof_map.C:1349
dof_id_type first_dof() const
Definition: dof_map.h:538
unsigned int idx(const ElemType type, const unsigned int nx, const unsigned int i, const unsigned int j)
void add_ghosting_functor(GhostingFunctor &ghosting_functor)
Definition: mesh_base.h:764
uint8_t dof_id_type
Definition: id_types.h:64
processor_id_type processor_id() const
dof_id_type _n_dfs
Definition: dof_map.h:1533
processor_id_type processor_id() const
Definition: dof_object.h:686
NodeConstraints _node_constraints
Definition: dof_map.h:1582
Defines the coupling between variables of a System.
void allgather(const T &send, std::vector< T > &recv) const
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2005