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
621  (static_cast<unsigned int>
622  (base_fe_type.order.get_order()),
623  FEInterface::max_order(base_fe_type,type),
624  "ERROR: Finite element "
625  << Utility::enum_to_string(base_fe_type.family)
626  << " on geometric element "
627  << Utility::enum_to_string(type)
628  << "\nonly supports FEInterface::max_order = "
629  << FEInterface::max_order(base_fe_type,type)
630  << ", not fe_type.order = "
631  << base_fe_type.order);
632 
633 # ifdef DEBUG
635  << "WARNING: Finite element "
636  << Utility::enum_to_string(base_fe_type.family)
637  << " on geometric element "
638  << Utility::enum_to_string(type) << std::endl
639  << "could not be p refined past FEInterface::max_order = "
640  << FEInterface::max_order(base_fe_type,type)
641  << std::endl;
642 # endif
643  elem->set_p_level(FEInterface::max_order(base_fe_type,type)
644  - base_fe_type.order);
645  }
646 #endif
647 
648  fe_type.order = static_cast<Order>(fe_type.order +
649  elem->p_level());
650 
651  // Allocate the vertex DOFs
652  for (unsigned int n=0; n<elem->n_nodes(); n++)
653  {
654  Node & node = elem->node_ref(n);
655 
656  if (elem->is_vertex(n))
657  {
658  const unsigned int old_node_dofs =
659  node.n_comp_group(sys_num, vg);
660 
661  const unsigned int vertex_dofs =
663  type, n),
664  old_node_dofs);
665 
666  // Some discontinuous FEs have no vertex dofs
667  if (vertex_dofs > old_node_dofs)
668  {
669  node.set_n_comp_group(sys_num, vg,
670  vertex_dofs);
671 
672  // Abusing dof_number to set a "this is a
673  // vertex" flag
674  node.set_vg_dof_base(sys_num, vg,
675  vertex_dofs);
676 
677  // libMesh::out << "sys_num,vg,old_node_dofs,vertex_dofs="
678  // << sys_num << ","
679  // << vg << ","
680  // << old_node_dofs << ","
681  // << vertex_dofs << '\n',
682  // node.debug_buffer();
683 
684  // libmesh_assert_equal_to (vertex_dofs, node.n_comp(sys_num, vg));
685  // libmesh_assert_equal_to (vertex_dofs, node.vg_dof_base(sys_num, vg));
686  }
687  }
688  }
689  } // done counting vertex dofs
690 
691  // count edge & face dofs next
692  elem_it = mesh.active_elements_begin();
693 
694  for ( ; elem_it != elem_end; ++elem_it)
695  {
696  Elem * elem = *elem_it;
697  libmesh_assert(elem);
698 
699  // Skip the numbering if this variable is
700  // not active on this element's subdomain
701  if (!vg_description.active_on_subdomain(elem->subdomain_id()))
702  continue;
703 
704  const ElemType type = elem->type();
705  const unsigned int dim = elem->dim();
706 
707  FEType fe_type = base_fe_type;
708  fe_type.order = static_cast<Order>(fe_type.order +
709  elem->p_level());
710 
711  // Allocate the edge and face DOFs
712  for (unsigned int n=0; n<elem->n_nodes(); n++)
713  {
714  Node & node = elem->node_ref(n);
715 
716  const unsigned int old_node_dofs =
717  node.n_comp_group(sys_num, vg);
718 
719  const unsigned int vertex_dofs = old_node_dofs?
720  cast_int<unsigned int>(node.vg_dof_base (sys_num,vg)):0;
721 
722  const unsigned int new_node_dofs =
723  FEInterface::n_dofs_at_node(dim, fe_type, type, n);
724 
725  // We've already allocated vertex DOFs
726  if (elem->is_vertex(n))
727  {
728  libmesh_assert_greater_equal (old_node_dofs, vertex_dofs);
729  // //if (vertex_dofs < new_node_dofs)
730  // libMesh::out << "sys_num,vg,old_node_dofs,vertex_dofs,new_node_dofs="
731  // << sys_num << ","
732  // << vg << ","
733  // << old_node_dofs << ","
734  // << vertex_dofs << ","
735  // << new_node_dofs << '\n',
736  // node.debug_buffer();
737 
738  libmesh_assert_greater_equal (vertex_dofs, new_node_dofs);
739  }
740  // We need to allocate the rest
741  else
742  {
743  // If this has no dofs yet, it needs no vertex
744  // dofs, so we just give it edge or face dofs
745  if (!old_node_dofs)
746  {
747  node.set_n_comp_group(sys_num, vg,
748  new_node_dofs);
749  // Abusing dof_number to set a "this has no
750  // vertex dofs" flag
751  if (new_node_dofs)
752  node.set_vg_dof_base(sys_num, vg, 0);
753  }
754 
755  // If this has dofs, but has no vertex dofs,
756  // it may still need more edge or face dofs if
757  // we're p-refined.
758  else if (vertex_dofs == 0)
759  {
760  if (new_node_dofs > old_node_dofs)
761  {
762  node.set_n_comp_group(sys_num, vg,
763  new_node_dofs);
764 
765  node.set_vg_dof_base(sys_num, vg,
766  vertex_dofs);
767  }
768  }
769  // If this is another element's vertex,
770  // add more (non-overlapping) edge/face dofs if
771  // necessary
772  else if (extra_hanging_dofs)
773  {
774  if (new_node_dofs > old_node_dofs - vertex_dofs)
775  {
776  node.set_n_comp_group(sys_num, vg,
777  vertex_dofs + new_node_dofs);
778 
779  node.set_vg_dof_base(sys_num, vg,
780  vertex_dofs);
781  }
782  }
783  // If this is another element's vertex, add any
784  // (overlapping) edge/face dofs if necessary
785  else
786  {
787  libmesh_assert_greater_equal (old_node_dofs, vertex_dofs);
788  if (new_node_dofs > old_node_dofs)
789  {
790  node.set_n_comp_group(sys_num, vg,
791  new_node_dofs);
792 
793  node.set_vg_dof_base (sys_num, vg,
794  vertex_dofs);
795  }
796  }
797  }
798  }
799  // Allocate the element DOFs
800  const unsigned int dofs_per_elem =
801  FEInterface::n_dofs_per_elem(dim, fe_type,
802  type);
803 
804  elem->set_n_comp_group(sys_num, vg, dofs_per_elem);
805 
806  }
807  } // end loop over variable groups
808 
809  // Calling DofMap::reinit() by itself makes little sense,
810  // so we won't bother with nonlocal DofObjects.
811  // Those will be fixed by distribute_dofs
812 
813  //------------------------------------------------------------
814  // Finally, clear all the current DOF indices
815  // (distribute_dofs expects them cleared!)
816  this->invalidate_dofs(mesh);
817 }
818 
819 
820 
822 {
823  const unsigned int sys_num = this->sys_number();
824 
825  // All the nodes
826  MeshBase::node_iterator node_it = mesh.nodes_begin();
827  const MeshBase::node_iterator node_end = mesh.nodes_end();
828 
829  for ( ; node_it != node_end; ++node_it)
830  (*node_it)->invalidate_dofs(sys_num);
831 
832  // All the elements
834  const MeshBase::element_iterator elem_end = mesh.active_elements_end();
835 
836  for ( ; elem_it != elem_end; ++elem_it)
837  (*elem_it)->invalidate_dofs(sys_num);
838 }
839 
840 
841 
843 {
844  // we don't want to clear
845  // the coupling matrix!
846  // It should not change...
847  //_dof_coupling->clear();
848  //
849  // But it would be inconsistent to leave our coupling settings
850  // through a clear()...
851  _dof_coupling = NULL;
852 
853  // Reset ghosting functor statuses
854  {
855  std::set<GhostingFunctor *>::iterator gf_it = this->coupling_functors_begin();
856  const std::set<GhostingFunctor *>::iterator gf_end = this->coupling_functors_end();
857  for (; gf_it != gf_end; ++gf_it)
858  {
859  GhostingFunctor *gf = *gf_it;
860  libmesh_assert(gf);
862  }
863  this->_coupling_functors.clear();
864 
865  // Go back to default coupling
866 
867  _default_coupling->set_dof_coupling(this->_dof_coupling);
868  _default_coupling->set_n_levels
869  (this->use_coupled_neighbor_dofs(this->_mesh));
870 
872  }
873 
874 
875  {
876  std::set<GhostingFunctor *>::iterator gf_it = this->algebraic_ghosting_functors_begin();
877  const std::set<GhostingFunctor *>::iterator gf_end = this->algebraic_ghosting_functors_end();
878  for (; gf_it != gf_end; ++gf_it)
879  {
880  GhostingFunctor *gf = *gf_it;
881  libmesh_assert(gf);
883  }
884  this->_algebraic_ghosting_functors.clear();
885 
886  // Go back to default send_list generation
887 
888  // _default_evaluating->set_dof_coupling(this->_dof_coupling);
889  _default_evaluating->set_n_levels(1);
891  }
892 
893  _variables.clear();
894  _variable_groups.clear();
895  _first_df.clear();
896  _end_df.clear();
897  _first_scalar_df.clear();
898  _send_list.clear();
899  this->clear_sparsity();
901 
902 #ifdef LIBMESH_ENABLE_AMR
903 
904  _dof_constraints.clear();
905  _stashed_dof_constraints.clear();
908  _n_old_dfs = 0;
909  _first_old_df.clear();
910  _end_old_df.clear();
911  _first_old_scalar_df.clear();
912 
913 #endif
914 
915  _matrices.clear();
916 
917  _n_dfs = 0;
918 }
919 
920 
921 
923 {
924  // This function must be run on all processors at once
925  parallel_object_only();
926 
927  // Log how long it takes to distribute the degrees of freedom
928  LOG_SCOPE("distribute_dofs()", "DofMap");
929 
930  libmesh_assert (mesh.is_prepared());
931 
932  const processor_id_type proc_id = this->processor_id();
933  const processor_id_type n_proc = this->n_processors();
934 
935  // libmesh_assert_greater (this->n_variables(), 0);
936  libmesh_assert_less (proc_id, n_proc);
937 
938  // re-init in case the mesh has changed
939  this->reinit(mesh);
940 
941  // By default distribute variables in a
942  // var-major fashion, but allow run-time
943  // specification
944  bool node_major_dofs = libMesh::on_command_line ("--node_major_dofs");
945 
946  // The DOF counter, will be incremented as we encounter
947  // new degrees of freedom
948  dof_id_type next_free_dof = 0;
949 
950  // Clear the send list before we rebuild it
951  _send_list.clear();
952 
953  // Set temporary DOF indices on this processor
954  if (node_major_dofs)
955  this->distribute_local_dofs_node_major (next_free_dof, mesh);
956  else
957  this->distribute_local_dofs_var_major (next_free_dof, mesh);
958 
959  // Get DOF counts on all processors
960  std::vector<dof_id_type> dofs_on_proc(n_proc, 0);
961  this->comm().allgather(next_free_dof, dofs_on_proc);
962 
963  // Resize and fill the _first_df and _end_df arrays
964 #ifdef LIBMESH_ENABLE_AMR
967 #endif
968 
969  _first_df.resize(n_proc);
970  _end_df.resize (n_proc);
971 
972  // Get DOF offsets
973  _first_df[0] = 0;
974  for (processor_id_type i=1; i < n_proc; ++i)
975  _first_df[i] = _end_df[i-1] = _first_df[i-1] + dofs_on_proc[i-1];
976  _end_df[n_proc-1] = _first_df[n_proc-1] + dofs_on_proc[n_proc-1];
977 
978  // Clear all the current DOF indices
979  // (distribute_dofs expects them cleared!)
980  this->invalidate_dofs(mesh);
981 
982  next_free_dof = _first_df[proc_id];
983 
984  // Set permanent DOF indices on this processor
985  if (node_major_dofs)
986  this->distribute_local_dofs_node_major (next_free_dof, mesh);
987  else
988  this->distribute_local_dofs_var_major (next_free_dof, mesh);
989 
990  libmesh_assert_equal_to (next_free_dof, _end_df[proc_id]);
991 
992  //------------------------------------------------------------
993  // At this point, all n_comp and dof_number values on local
994  // DofObjects should be correct, but a DistributedMesh might have
995  // incorrect values on non-local DofObjects. Let's request the
996  // correct values from each other processor.
997 
998  if (this->n_processors() > 1)
999  {
1001  mesh.nodes_end(),
1003 
1005  mesh.elements_end(),
1007  }
1008 
1009 #ifdef DEBUG
1010  {
1011  const unsigned int
1012  sys_num = this->sys_number();
1013 
1014  // Processors should all agree on DoF ids for the newly numbered
1015  // system.
1017 
1018  // DoF processor ids should match DofObject processor ids
1019  MeshBase::const_node_iterator node_it = mesh.nodes_begin();
1020  const MeshBase::const_node_iterator node_end = mesh.nodes_end();
1021  for ( ; node_it != node_end; ++node_it)
1022  {
1023  DofObject const * const dofobj = *node_it;
1024  const processor_id_type obj_proc_id = dofobj->processor_id();
1025 
1026  for (unsigned int v=0; v != dofobj->n_vars(sys_num); ++v)
1027  for (unsigned int c=0; c != dofobj->n_comp(sys_num,v); ++c)
1028  {
1029  const dof_id_type dofid = dofobj->dof_number(sys_num,v,c);
1030  libmesh_assert_greater_equal (dofid, this->first_dof(obj_proc_id));
1031  libmesh_assert_less (dofid, this->end_dof(obj_proc_id));
1032  }
1033  }
1034 
1036  const MeshBase::const_element_iterator elem_end = mesh.elements_end();
1037  for ( ; elem_it != elem_end; ++elem_it)
1038  {
1039  DofObject const * const dofobj = *elem_it;
1040  const processor_id_type obj_proc_id = dofobj->processor_id();
1041 
1042  for (unsigned int v=0; v != dofobj->n_vars(sys_num); ++v)
1043  for (unsigned int c=0; c != dofobj->n_comp(sys_num,v); ++c)
1044  {
1045  const dof_id_type dofid = dofobj->dof_number(sys_num,v,c);
1046  libmesh_assert_greater_equal (dofid, this->first_dof(obj_proc_id));
1047  libmesh_assert_less (dofid, this->end_dof(obj_proc_id));
1048  }
1049  }
1050  }
1051 #endif
1052 
1053  // Set the total number of degrees of freedom, then start finding
1054  // SCALAR degrees of freedom
1055 #ifdef LIBMESH_ENABLE_AMR
1056  _n_old_dfs = _n_dfs;
1058 #endif
1059  _n_dfs = _end_df[n_proc-1];
1060  _first_scalar_df.clear();
1062  dof_id_type current_SCALAR_dof_index = n_dofs() - n_SCALAR_dofs();
1063 
1064  // Calculate and cache the initial DoF indices for SCALAR variables.
1065  // This is an O(N_vars) calculation so we want to do it once per
1066  // renumbering rather than once per SCALAR_dof_indices() call
1067 
1068  for (unsigned int v=0; v<this->n_variables(); v++)
1069  if(this->variable(v).type().family == SCALAR)
1070  {
1071  _first_scalar_df[v] = current_SCALAR_dof_index;
1072  current_SCALAR_dof_index += this->variable(v).type().order.get_order();
1073  }
1074 
1075  // Allow our GhostingFunctor objects to reinit if necessary
1076  {
1077  std::set<GhostingFunctor *>::iterator gf_it = this->algebraic_ghosting_functors_begin();
1078  const std::set<GhostingFunctor *>::iterator gf_end = this->algebraic_ghosting_functors_end();
1079  for (; gf_it != gf_end; ++gf_it)
1080  {
1081  GhostingFunctor *gf = *gf_it;
1082  libmesh_assert(gf);
1083  gf->dofmap_reinit();
1084  }
1085  }
1086 
1087  {
1088  std::set<GhostingFunctor *>::iterator gf_it = this->coupling_functors_begin();
1089  const std::set<GhostingFunctor *>::iterator gf_end = this->coupling_functors_end();
1090  for (; gf_it != gf_end; ++gf_it)
1091  {
1092  GhostingFunctor *gf = *gf_it;
1093  libmesh_assert(gf);
1094  gf->dofmap_reinit();
1095  }
1096  }
1097 
1098  // Note that in the add_neighbors_to_send_list nodes on processor
1099  // boundaries that are shared by multiple elements are added for
1100  // each element.
1101  this->add_neighbors_to_send_list(mesh);
1102 
1103  // Here we used to clean up that data structure; now System and
1104  // EquationSystems call that for us, after we've added constraint
1105  // dependencies to the send_list too.
1106  // this->sort_send_list ();
1107 }
1108 
1109 
1110 void DofMap::local_variable_indices(std::vector<dof_id_type> & idx,
1111  const MeshBase & mesh,
1112  unsigned int var_num) const
1113 {
1114  const unsigned int sys_num = this->sys_number();
1115 
1116  // If this isn't a SCALAR variable, we need to find all its field
1117  // dofs on the mesh
1118  if (this->variable_type(var_num).family != SCALAR)
1119  {
1122 
1123  for ( ; elem_it != elem_end; ++elem_it)
1124  {
1125  // Only count dofs connected to active
1126  // elements on this processor.
1127  Elem * elem = *elem_it;
1128  const unsigned int n_nodes = elem->n_nodes();
1129 
1130  // First get any new nodal DOFS
1131  for (unsigned int n=0; n<n_nodes; n++)
1132  {
1133  Node & node = elem->node_ref(n);
1134 
1135  if (node.processor_id() < this->processor_id())
1136  continue;
1137 
1138  const unsigned int n_comp = node.n_comp(sys_num, var_num);
1139  for(unsigned int i=0; i<n_comp; i++)
1140  {
1141  const dof_id_type index = node.dof_number(sys_num,var_num,i);
1142  libmesh_assert_greater_equal (index, this->first_dof());
1143  libmesh_assert_less (index, this->end_dof());
1144 
1145  if (idx.empty() || index > idx.back())
1146  idx.push_back(index);
1147  }
1148  }
1149 
1150  // Next get any new element DOFS
1151  const unsigned int n_comp = elem->n_comp(sys_num, var_num);
1152  for(unsigned int i=0; i<n_comp; i++)
1153  {
1154  const dof_id_type index = elem->dof_number(sys_num,var_num,i);
1155  if (idx.empty() || index > idx.back())
1156  idx.push_back(index);
1157  }
1158  } // done looping over elements
1159 
1160 
1161  // we may have missed assigning DOFs to nodes that we own
1162  // but to which we have no connected elements matching our
1163  // variable restriction criterion. this will happen, for example,
1164  // if variable V is restricted to subdomain S. We may not own
1165  // any elements which live in S, but we may own nodes which are
1166  // *connected* to elements which do. in this scenario these nodes
1167  // will presently have unnumbered DOFs. we need to take care of
1168  // them here since we own them and no other processor will touch them.
1169  {
1171  const MeshBase::const_node_iterator node_end = mesh.local_nodes_end();
1172 
1173  for (; node_it != node_end; ++node_it)
1174  {
1175  Node * node = *node_it;
1176  libmesh_assert(node);
1177 
1178  const unsigned int n_comp = node->n_comp(sys_num, var_num);
1179  for(unsigned int i=0; i<n_comp; i++)
1180  {
1181  const dof_id_type index = node->dof_number(sys_num,var_num,i);
1182  if (idx.empty() || index > idx.back())
1183  idx.push_back(index);
1184  }
1185  }
1186  }
1187  }
1188  // Otherwise, count up the SCALAR dofs, if we're on the processor
1189  // that holds this SCALAR variable
1190  else if ( this->processor_id() == (this->n_processors()-1) )
1191  {
1192  std::vector<dof_id_type> di_scalar;
1193  this->SCALAR_dof_indices(di_scalar,var_num);
1194  idx.insert( idx.end(), di_scalar.begin(), di_scalar.end());
1195  }
1196 }
1197 
1198 
1200  MeshBase & mesh)
1201 {
1202  const unsigned int sys_num = this->sys_number();
1203  const unsigned int n_var_groups = this->n_variable_groups();
1204 
1205  //-------------------------------------------------------------------------
1206  // First count and assign temporary numbers to local dofs
1208  const MeshBase::element_iterator elem_end = mesh.active_local_elements_end();
1209 
1210  for ( ; elem_it != elem_end; ++elem_it)
1211  {
1212  // Only number dofs connected to active
1213  // elements on this processor.
1214  Elem * elem = *elem_it;
1215  const unsigned int n_nodes = elem->n_nodes();
1216 
1217  // First number the nodal DOFS
1218  for (unsigned int n=0; n<n_nodes; n++)
1219  {
1220  Node & node = elem->node_ref(n);
1221 
1222  for (unsigned vg=0; vg<n_var_groups; vg++)
1223  {
1224  const VariableGroup & vg_description(this->variable_group(vg));
1225 
1226  if( (vg_description.type().family != SCALAR) &&
1227  (vg_description.active_on_subdomain(elem->subdomain_id())) )
1228  {
1229  // assign dof numbers (all at once) if this is
1230  // our node and if they aren't already there
1231  if ((node.n_comp_group(sys_num,vg) > 0) &&
1232  (node.processor_id() == this->processor_id()) &&
1233  (node.vg_dof_base(sys_num,vg) ==
1235  {
1236  node.set_vg_dof_base(sys_num, vg,
1237  next_free_dof);
1238  next_free_dof += (vg_description.n_variables()*
1239  node.n_comp_group(sys_num,vg));
1240  //node.debug_buffer();
1241  }
1242  }
1243  }
1244  }
1245 
1246  // Now number the element DOFS
1247  for (unsigned vg=0; vg<n_var_groups; vg++)
1248  {
1249  const VariableGroup & vg_description(this->variable_group(vg));
1250 
1251  if ( (vg_description.type().family != SCALAR) &&
1252  (vg_description.active_on_subdomain(elem->subdomain_id())) )
1253  if (elem->n_comp_group(sys_num,vg) > 0)
1254  {
1255  libmesh_assert_equal_to (elem->vg_dof_base(sys_num,vg),
1257 
1258  elem->set_vg_dof_base(sys_num,
1259  vg,
1260  next_free_dof);
1261 
1262  next_free_dof += (vg_description.n_variables()*
1263  elem->n_comp(sys_num,vg));
1264  }
1265  }
1266  } // done looping over elements
1267 
1268 
1269  // we may have missed assigning DOFs to nodes that we own
1270  // but to which we have no connected elements matching our
1271  // variable restriction criterion. this will happen, for example,
1272  // if variable V is restricted to subdomain S. We may not own
1273  // any elements which live in S, but we may own nodes which are
1274  // *connected* to elements which do. in this scenario these nodes
1275  // will presently have unnumbered DOFs. we need to take care of
1276  // them here since we own them and no other processor will touch them.
1277  {
1278  MeshBase::node_iterator node_it = mesh.local_nodes_begin();
1279  const MeshBase::node_iterator node_end = mesh.local_nodes_end();
1280 
1281  for (; node_it != node_end; ++node_it)
1282  {
1283  Node * node = *node_it;
1284  libmesh_assert(node);
1285 
1286  for (unsigned vg=0; vg<n_var_groups; vg++)
1287  {
1288  const VariableGroup & vg_description(this->variable_group(vg));
1289 
1290  if (node->n_comp_group(sys_num,vg))
1291  if (node->vg_dof_base(sys_num,vg) == DofObject::invalid_id)
1292  {
1293  node->set_vg_dof_base (sys_num,
1294  vg,
1295  next_free_dof);
1296 
1297  next_free_dof += (vg_description.n_variables()*
1298  node->n_comp(sys_num,vg));
1299  }
1300  }
1301  }
1302  }
1303 
1304  // Finally, count up the SCALAR dofs
1305  this->_n_SCALAR_dofs = 0;
1306  for (unsigned vg=0; vg<n_var_groups; vg++)
1307  {
1308  const VariableGroup & vg_description(this->variable_group(vg));
1309 
1310  if( vg_description.type().family == SCALAR )
1311  {
1312  this->_n_SCALAR_dofs += (vg_description.n_variables()*
1313  vg_description.type().order.get_order());
1314  continue;
1315  }
1316  }
1317 
1318  // Only increment next_free_dof if we're on the processor
1319  // that holds this SCALAR variable
1320  if ( this->processor_id() == (this->n_processors()-1) )
1321  next_free_dof += _n_SCALAR_dofs;
1322 
1323 #ifdef DEBUG
1324  {
1325  // libMesh::out << "next_free_dof=" << next_free_dof << std::endl
1326  // << "_n_SCALAR_dofs=" << _n_SCALAR_dofs << std::endl;
1327 
1328  // Make sure we didn't miss any nodes
1329  MeshTools::libmesh_assert_valid_procids<Node>(mesh);
1330 
1331  MeshBase::node_iterator node_it = mesh.local_nodes_begin();
1332  const MeshBase::node_iterator node_end = mesh.local_nodes_end();
1333  for (; node_it != node_end; ++node_it)
1334  {
1335  Node * obj = *node_it;
1336  libmesh_assert(obj);
1337  unsigned int n_var_g = obj->n_var_groups(this->sys_number());
1338  for (unsigned int vg=0; vg != n_var_g; ++vg)
1339  {
1340  unsigned int n_comp_g =
1341  obj->n_comp_group(this->sys_number(), vg);
1342  dof_id_type my_first_dof = n_comp_g ?
1343  obj->vg_dof_base(this->sys_number(), vg) : 0;
1344  libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
1345  }
1346  }
1347  }
1348 #endif // DEBUG
1349 }
1350 
1351 
1352 
1354  MeshBase & mesh)
1355 {
1356  const unsigned int sys_num = this->sys_number();
1357  const unsigned int n_var_groups = this->n_variable_groups();
1358 
1359  //-------------------------------------------------------------------------
1360  // First count and assign temporary numbers to local dofs
1361  for (unsigned vg=0; vg<n_var_groups; vg++)
1362  {
1363  const VariableGroup & vg_description(this->variable_group(vg));
1364 
1365  const unsigned int n_vars_in_group = vg_description.n_variables();
1366 
1367  // Skip the SCALAR dofs
1368  if (vg_description.type().family == SCALAR)
1369  continue;
1370 
1372  const MeshBase::element_iterator elem_end = mesh.active_local_elements_end();
1373 
1374  for ( ; elem_it != elem_end; ++elem_it)
1375  {
1376  // Only number dofs connected to active
1377  // elements on this processor.
1378  Elem * elem = *elem_it;
1379 
1380  // ... and only variables which are active on
1381  // on this element's subdomain
1382  if (!vg_description.active_on_subdomain(elem->subdomain_id()))
1383  continue;
1384 
1385  const unsigned int n_nodes = elem->n_nodes();
1386 
1387  // First number the nodal DOFS
1388  for (unsigned int n=0; n<n_nodes; n++)
1389  {
1390  Node & node = elem->node_ref(n);
1391 
1392  // assign dof numbers (all at once) if this is
1393  // our node and if they aren't already there
1394  if ((node.n_comp_group(sys_num,vg) > 0) &&
1395  (node.processor_id() == this->processor_id()) &&
1396  (node.vg_dof_base(sys_num,vg) ==
1398  {
1399  node.set_vg_dof_base(sys_num, vg, next_free_dof);
1400 
1401  next_free_dof += (n_vars_in_group*
1402  node.n_comp_group(sys_num,vg));
1403  }
1404  }
1405 
1406  // Now number the element DOFS
1407  if (elem->n_comp_group(sys_num,vg) > 0)
1408  {
1409  libmesh_assert_equal_to (elem->vg_dof_base(sys_num,vg),
1411 
1412  elem->set_vg_dof_base(sys_num,
1413  vg,
1414  next_free_dof);
1415 
1416  next_free_dof += (n_vars_in_group*
1417  elem->n_comp_group(sys_num,vg));
1418  }
1419  } // end loop on elements
1420 
1421  // we may have missed assigning DOFs to nodes that we own
1422  // but to which we have no connected elements matching our
1423  // variable restriction criterion. this will happen, for example,
1424  // if variable V is restricted to subdomain S. We may not own
1425  // any elements which live in S, but we may own nodes which are
1426  // *connected* to elements which do. in this scenario these nodes
1427  // will presently have unnumbered DOFs. we need to take care of
1428  // them here since we own them and no other processor will touch them.
1429  {
1430  MeshBase::node_iterator node_it = mesh.local_nodes_begin();
1431  const MeshBase::node_iterator node_end = mesh.local_nodes_end();
1432 
1433  for (; node_it != node_end; ++node_it)
1434  {
1435  Node * node = *node_it;
1436  libmesh_assert(node);
1437 
1438  if (node->n_comp_group(sys_num,vg))
1439  if (node->vg_dof_base(sys_num,vg) == DofObject::invalid_id)
1440  {
1441  node->set_vg_dof_base (sys_num,
1442  vg,
1443  next_free_dof);
1444 
1445  next_free_dof += (n_vars_in_group*
1446  node->n_comp_group(sys_num,vg));
1447  }
1448  }
1449  }
1450  } // end loop on variable groups
1451 
1452  // Finally, count up the SCALAR dofs
1453  this->_n_SCALAR_dofs = 0;
1454  for (unsigned vg=0; vg<n_var_groups; vg++)
1455  {
1456  const VariableGroup & vg_description(this->variable_group(vg));
1457 
1458  if( vg_description.type().family == SCALAR )
1459  {
1460  this->_n_SCALAR_dofs += (vg_description.n_variables()*
1461  vg_description.type().order.get_order());
1462  continue;
1463  }
1464  }
1465 
1466  // Only increment next_free_dof if we're on the processor
1467  // that holds this SCALAR variable
1468  if ( this->processor_id() == (this->n_processors()-1) )
1469  next_free_dof += _n_SCALAR_dofs;
1470 
1471 #ifdef DEBUG
1472  {
1473  // Make sure we didn't miss any nodes
1474  MeshTools::libmesh_assert_valid_procids<Node>(mesh);
1475 
1476  MeshBase::node_iterator node_it = mesh.local_nodes_begin();
1477  const MeshBase::node_iterator node_end = mesh.local_nodes_end();
1478  for (; node_it != node_end; ++node_it)
1479  {
1480  Node * obj = *node_it;
1481  libmesh_assert(obj);
1482  unsigned int n_var_g = obj->n_var_groups(this->sys_number());
1483  for (unsigned int vg=0; vg != n_var_g; ++vg)
1484  {
1485  unsigned int n_comp_g =
1486  obj->n_comp_group(this->sys_number(), vg);
1487  dof_id_type my_first_dof = n_comp_g ?
1488  obj->vg_dof_base(this->sys_number(), vg) : 0;
1489  libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
1490  }
1491  }
1492  }
1493 #endif // DEBUG
1494 }
1495 
1496 
1497 
1499  (GhostingFunctor::map_type & elements_to_ghost,
1500  std::set<CouplingMatrix*> & temporary_coupling_matrices,
1501  const std::set<GhostingFunctor *>::iterator & gf_begin,
1502  const std::set<GhostingFunctor *>::iterator & gf_end,
1503  const MeshBase::const_element_iterator & elems_begin,
1504  const MeshBase::const_element_iterator & elems_end,
1506 {
1507  std::set<GhostingFunctor *>::iterator gf_it = gf_begin;
1508  for (; gf_it != gf_end; ++gf_it)
1509  {
1510  GhostingFunctor::map_type more_elements_to_ghost;
1511 
1512  GhostingFunctor *gf = *gf_it;
1513  libmesh_assert(gf);
1514  (*gf)(elems_begin, elems_end, p, more_elements_to_ghost);
1515 
1516  GhostingFunctor::map_type::iterator metg_it = more_elements_to_ghost.begin();
1517  const GhostingFunctor::map_type::iterator metg_end = more_elements_to_ghost.end();
1518  for (; metg_it != metg_end; ++metg_it)
1519  {
1520  GhostingFunctor::map_type::iterator existing_it =
1521  elements_to_ghost.find (metg_it->first);
1522  if (existing_it == elements_to_ghost.end())
1523  elements_to_ghost.insert(*metg_it);
1524  else
1525  {
1526  if (existing_it->second)
1527  {
1528  if (metg_it->second)
1529  {
1530  // If this isn't already a temporary
1531  // then we need to make one so we'll
1532  // have a non-const matrix to merge
1533  if (temporary_coupling_matrices.empty() ||
1534  temporary_coupling_matrices.find
1535  (const_cast<CouplingMatrix*>(existing_it->second)) ==
1536  temporary_coupling_matrices.end())
1537  {
1538  CouplingMatrix *cm = new CouplingMatrix(*existing_it->second);
1539  temporary_coupling_matrices.insert(cm);
1540  existing_it->second = cm;
1541  }
1542  const_cast<CouplingMatrix&>(*existing_it->second) &= *metg_it->second;
1543  }
1544  else
1545  {
1546  // Any existing_it matrix merged with a full
1547  // matrix (symbolized as nullptr) gives another
1548  // full matrix (symbolizable as nullptr).
1549 
1550  // So if existing_it->second is a temporary then
1551  // we don't need it anymore; we might as well
1552  // remove it to keep the set of temporaries
1553  // small.
1554  std::set<CouplingMatrix*>::iterator temp_it =
1555  temporary_coupling_matrices.find
1556  (const_cast<CouplingMatrix*>(existing_it->second));
1557  if (temp_it !=
1558  temporary_coupling_matrices.end())
1559  temporary_coupling_matrices.erase(temp_it);
1560 
1561  existing_it->second = libmesh_nullptr;
1562  }
1563  }
1564  // else we have a nullptr already, then we have a full
1565  // coupling matrix, already, and merging with anything
1566  // else won't change that, so we're done.
1567  }
1568  }
1569  }
1570 }
1571 
1572 
1573 
1575 {
1576  LOG_SCOPE("add_neighbors_to_send_list()", "DofMap");
1577 
1578  const unsigned int n_var = this->n_variables();
1579 
1580  MeshBase::const_element_iterator local_elem_it
1581  = mesh.active_local_elements_begin();
1582  const MeshBase::const_element_iterator local_elem_end
1583  = mesh.active_local_elements_end();
1584 
1585  GhostingFunctor::map_type elements_to_send;
1586 
1587  // Man, I wish we had guaranteed unique_ptr availability...
1588  std::set<CouplingMatrix*> temporary_coupling_matrices;
1589 
1590  // We need to add dofs to the send list if they've been directly
1591  // requested by an algebraic ghosting functor or they've been
1592  // indirectly requested by a coupling functor.
1594  (elements_to_send,
1595  temporary_coupling_matrices,
1598  local_elem_it, local_elem_end, mesh.processor_id());
1599 
1601  (elements_to_send,
1602  temporary_coupling_matrices,
1603  this->coupling_functors_begin(),
1604  this->coupling_functors_end(),
1605  local_elem_it, local_elem_end, mesh.processor_id());
1606 
1607  // Making a list of non-zero coupling matrix columns is an
1608  // O(N_var^2) operation. We cache it so we only have to do it once
1609  // per CouplingMatrix and not once per element.
1610  std::map<const CouplingMatrix*, std::vector<unsigned int> >
1611  column_variable_lists;
1612 
1613  GhostingFunctor::map_type::iterator etg_it = elements_to_send.begin();
1614  const GhostingFunctor::map_type::iterator etg_end = elements_to_send.end();
1615  for (; etg_it != etg_end; ++etg_it)
1616  {
1617  const Elem * const partner = etg_it->first;
1618 
1619  // We asked ghosting functors not to give us local elements
1620  libmesh_assert_not_equal_to
1621  (partner->processor_id(), this->processor_id());
1622 
1623  const CouplingMatrix *ghost_coupling = etg_it->second;
1624 
1625  // Loop over any present coupling matrix column variables if we
1626  // have a coupling matrix, or just add all variables to
1627  // send_list if not.
1628  if (ghost_coupling)
1629  {
1630  libmesh_assert_equal_to (ghost_coupling->size(), n_var);
1631 
1632  // Try to find a cached list of column variables.
1633  std::map<const CouplingMatrix*,
1634  std::vector<unsigned int> >::const_iterator
1635  column_variable_list =
1636  column_variable_lists.find(ghost_coupling);
1637 
1638  // If we didn't find it, then we need to create it.
1639  if (column_variable_list == column_variable_lists.end())
1640  {
1641  std::pair<
1642  std::map<const CouplingMatrix*,
1643  std::vector<unsigned int> >::iterator,
1644  bool>
1645  inserted_variable_list_pair =
1646  column_variable_lists.insert
1647  (std::pair<const CouplingMatrix*,
1648  std::vector<unsigned int> >
1649  (ghost_coupling, std::vector<unsigned int>()));
1650  column_variable_list = inserted_variable_list_pair.first;
1651 
1652  std::vector<unsigned int> & new_variable_list =
1653  inserted_variable_list_pair.first->second;
1654 
1655  std::vector<unsigned char> has_variable(n_var, false);
1656 
1657  for (unsigned int vi = 0; vi != n_var; ++vi)
1658  {
1659  ConstCouplingRow ccr(vi, *ghost_coupling);
1660 
1661  for (ConstCouplingRow::const_iterator it = ccr.begin(),
1662  end = ccr.end();
1663  it != end; ++it)
1664  {
1665  const unsigned int vj = *it;
1666  has_variable[vj] = true;
1667  }
1668  }
1669  for (unsigned int vj = 0; vj != n_var; ++vj)
1670  {
1671  if (has_variable[vj])
1672  new_variable_list.push_back(vj);
1673  }
1674  }
1675 
1676  const std::vector<unsigned int> & variable_list =
1677  column_variable_list->second;
1678 
1679  for (std::size_t j=0; j != variable_list.size(); ++j)
1680  {
1681  const unsigned int vj=variable_list[j];
1682 
1683  std::vector<dof_id_type> di;
1684  this->dof_indices (partner, di, vj);
1685 
1686  // Insert the remote DOF indices into the send list
1687  for (std::size_t j=0; j != di.size(); ++j)
1688  if (di[j] < this->first_dof() ||
1689  di[j] >= this->end_dof())
1690  _send_list.push_back(di[j]);
1691  }
1692  }
1693  else
1694  {
1695  std::vector<dof_id_type> di;
1696  this->dof_indices (partner, di);
1697 
1698  // Insert the remote DOF indices into the send list
1699  for (std::size_t j=0; j != di.size(); ++j)
1700  if (di[j] < this->first_dof() ||
1701  di[j] >= this->end_dof())
1702  _send_list.push_back(di[j]);
1703  }
1704 
1705  }
1706 
1707  // We're now done with any merged coupling matrices we had to create.
1708  for (std::set<CouplingMatrix*>::iterator
1709  it = temporary_coupling_matrices.begin(),
1710  end = temporary_coupling_matrices.begin();
1711  it != end; ++it)
1712  delete *it;
1713 
1714  //-------------------------------------------------------------------------
1715  // Our coupling functors added dofs from neighboring elements to the
1716  // send list, but we may still need to add non-local dofs from local
1717  // elements.
1718  //-------------------------------------------------------------------------
1719 
1720  // Loop over the active local elements, adding all active elements
1721  // that neighbor an active local element to the send list.
1722  for ( ; local_elem_it != local_elem_end; ++local_elem_it)
1723  {
1724  const Elem * elem = *local_elem_it;
1725 
1726  std::vector<dof_id_type> di;
1727  this->dof_indices (elem, di);
1728 
1729  // Insert the remote DOF indices into the send list
1730  for (std::size_t j=0; j != di.size(); ++j)
1731  if (di[j] < this->first_dof() ||
1732  di[j] >= this->end_dof())
1733  _send_list.push_back(di[j]);
1734  }
1735 }
1736 
1737 
1738 
1740 {
1741  LOG_SCOPE("prepare_send_list()", "DofMap");
1742 
1743  // Check to see if we have any extra stuff to add to the send_list
1745  {
1746  if (_augment_send_list)
1747  {
1748  libmesh_here();
1749  libMesh::out << "WARNING: You have specified both an extra send list function and object.\n"
1750  << " Are you sure this is what you meant to do??"
1751  << std::endl;
1752  }
1753 
1755  }
1756 
1757  if (_augment_send_list)
1759 
1760  // First sort the send list. After this
1761  // duplicated elements will be adjacent in the
1762  // vector
1763  std::sort(_send_list.begin(), _send_list.end());
1764 
1765  // Now use std::unique to remove duplicate entries
1766  std::vector<dof_id_type>::iterator new_end =
1767  std::unique (_send_list.begin(), _send_list.end());
1768 
1769  // Remove the end of the send_list. Use the "swap trick"
1770  // from Effective STL
1771  std::vector<dof_id_type> (_send_list.begin(), new_end).swap (_send_list);
1772 }
1773 
1774 
1775 void DofMap::set_implicit_neighbor_dofs(bool implicit_neighbor_dofs)
1776 {
1778  _implicit_neighbor_dofs = implicit_neighbor_dofs;
1779 }
1780 
1781 
1783 {
1784  // If we were asked on the command line, then we need to
1785  // include sensitivities between neighbor degrees of freedom
1786  bool implicit_neighbor_dofs =
1787  libMesh::on_command_line ("--implicit_neighbor_dofs");
1788 
1789  // If the user specifies --implicit_neighbor_dofs 0, then
1790  // presumably he knows what he is doing and we won't try to
1791  // automatically turn it on even when all the variables are
1792  // discontinuous.
1793  if (implicit_neighbor_dofs)
1794  {
1795  // No flag provided defaults to 'true'
1796  int flag = 1;
1797  flag = libMesh::command_line_next ("--implicit_neighbor_dofs", flag);
1798 
1799  if (!flag)
1800  {
1801  // The user said --implicit_neighbor_dofs 0, so he knows
1802  // what he is doing and really doesn't want it.
1803  return false;
1804  }
1805  }
1806 
1807  // Possibly override the commandline option, if set_implicit_neighbor_dofs
1808  // has been called.
1810  {
1811  implicit_neighbor_dofs = _implicit_neighbor_dofs;
1812 
1813  // Again, if the user explicitly says implicit_neighbor_dofs = false,
1814  // then we return here.
1815  if (!implicit_neighbor_dofs)
1816  return false;
1817  }
1818 
1819  // Look at all the variables in this system. If every one is
1820  // discontinuous then the user must be doing DG/FVM, so be nice
1821  // and force implicit_neighbor_dofs=true.
1822  {
1823  bool all_discontinuous_dofs = true;
1824 
1825  for (unsigned int var=0; var<this->n_variables(); var++)
1826  if (FEAbstract::build (mesh.mesh_dimension(),
1827  this->variable_type(var))->get_continuity() != DISCONTINUOUS)
1828  all_discontinuous_dofs = false;
1829 
1830  if (all_discontinuous_dofs)
1831  implicit_neighbor_dofs = true;
1832  }
1833 
1834  return implicit_neighbor_dofs;
1835 }
1836 
1837 
1838 
1840 {
1841  _sp = this->build_sparsity(mesh);
1842 
1843  // It is possible that some \p SparseMatrix implementations want to
1844  // see it. Let them see it before we throw it away.
1845  std::vector<SparseMatrix<Number> *>::const_iterator
1846  pos = _matrices.begin(),
1847  end = _matrices.end();
1848 
1849  // If we need the full sparsity pattern, then we share a view of its
1850  // arrays, and we pass it in to the matrices.
1852  {
1853  _n_nz = &_sp->n_nz;
1854  _n_oz = &_sp->n_oz;
1855 
1856  for (; pos != end; ++pos)
1857  (*pos)->update_sparsity_pattern (_sp->sparsity_pattern);
1858  }
1859  // If we don't need the full sparsity pattern anymore, steal the
1860  // arrays we do need and free the rest of the memory
1861  else
1862  {
1863  if (!_n_nz)
1864  _n_nz = new std::vector<dof_id_type>();
1865  _n_nz->swap(_sp->n_nz);
1866  if (!_n_oz)
1867  _n_oz = new std::vector<dof_id_type>();
1868  _n_oz->swap(_sp->n_oz);
1869 
1870  _sp.reset();
1871  }
1872 }
1873 
1874 
1875 
1877 {
1879  {
1880  libmesh_assert(_sp.get());
1881  libmesh_assert(!_n_nz || _n_nz == &_sp->n_nz);
1882  libmesh_assert(!_n_oz || _n_oz == &_sp->n_oz);
1883  _sp.reset();
1884  }
1885  else
1886  {
1887  libmesh_assert(!_sp.get());
1888  delete _n_nz;
1889  delete _n_oz;
1890  }
1893 }
1894 
1895 
1896 
1898  (GhostingFunctor & coupling_functor)
1899 {
1900  _coupling_functors.insert(&coupling_functor);
1901  _mesh.add_ghosting_functor(coupling_functor);
1902 }
1903 
1904 
1905 
1907  (GhostingFunctor & coupling_functor)
1908 {
1909  libmesh_assert(_coupling_functors.count(&coupling_functor));
1910  _coupling_functors.erase(&coupling_functor);
1911  _mesh.remove_ghosting_functor(coupling_functor);
1912 }
1913 
1914 
1915 
1917  (GhostingFunctor & algebraic_ghosting_functor)
1918 {
1919  _algebraic_ghosting_functors.insert(&algebraic_ghosting_functor);
1920  _mesh.add_ghosting_functor(algebraic_ghosting_functor);
1921 }
1922 
1923 
1924 
1926  (GhostingFunctor & algebraic_ghosting_functor)
1927 {
1929  (&algebraic_ghosting_functor));
1930  _algebraic_ghosting_functors.erase(&algebraic_ghosting_functor);
1931  _mesh.remove_ghosting_functor(algebraic_ghosting_functor);
1932 }
1933 
1934 
1935 
1937  const std::vector<dof_id_type> & dof_indices_in,
1938  DenseVectorBase<Number> & Ue) const
1939 {
1940 #ifdef LIBMESH_ENABLE_AMR
1941 
1942  // Trivial mapping
1943  libmesh_assert_equal_to (dof_indices_in.size(), Ue.size());
1944  bool has_constrained_dofs = false;
1945 
1946  for (unsigned int il=0;
1947  il != cast_int<unsigned int>(dof_indices_in.size()); il++)
1948  {
1949  const dof_id_type ig = dof_indices_in[il];
1950 
1951  if (this->is_constrained_dof (ig)) has_constrained_dofs = true;
1952 
1953  libmesh_assert_less (ig, Ug.size());
1954 
1955  Ue.el(il) = Ug(ig);
1956  }
1957 
1958  // If the element has any constrained DOFs then we need
1959  // to account for them in the mapping. This will handle
1960  // the case that the input vector is not constrained.
1961  if (has_constrained_dofs)
1962  {
1963  // Copy the input DOF indices.
1964  std::vector<dof_id_type> constrained_dof_indices(dof_indices_in);
1965 
1968 
1969  this->build_constraint_matrix_and_vector (C, H, constrained_dof_indices);
1970 
1971  libmesh_assert_equal_to (dof_indices_in.size(), C.m());
1972  libmesh_assert_equal_to (constrained_dof_indices.size(), C.n());
1973 
1974  // zero-out Ue
1975  Ue.zero();
1976 
1977  // compute Ue = C Ug, with proper mapping.
1978  const unsigned int n_original_dofs =
1979  cast_int<unsigned int>(dof_indices_in.size());
1980  for (unsigned int i=0; i != n_original_dofs; i++)
1981  {
1982  Ue.el(i) = H(i);
1983 
1984  const unsigned int n_constrained =
1985  cast_int<unsigned int>(constrained_dof_indices.size());
1986  for (unsigned int j=0; j<n_constrained; j++)
1987  {
1988  const dof_id_type jg = constrained_dof_indices[j];
1989 
1990  // If Ug is a serial or ghosted vector, then this assert is
1991  // overzealous. If Ug is a parallel vector, then this assert
1992  // is redundant.
1993  // libmesh_assert ((jg >= Ug.first_local_index()) &&
1994  // (jg < Ug.last_local_index()));
1995 
1996  Ue.el(i) += C(i,j)*Ug(jg);
1997  }
1998  }
1999  }
2000 
2001 #else
2002 
2003  // Trivial mapping
2004 
2005  const unsigned int n_original_dofs =
2006  cast_int<unsigned int>(dof_indices_in.size());
2007 
2008  libmesh_assert_equal_to (n_original_dofs, Ue.size());
2009 
2010  for (unsigned int il=0; il<n_original_dofs; il++)
2011  {
2012  const dof_id_type ig = dof_indices_in[il];
2013 
2014  libmesh_assert ((ig >= Ug.first_local_index()) && (ig < Ug.last_local_index()));
2015 
2016  Ue.el(il) = Ug(ig);
2017  }
2018 
2019 #endif
2020 }
2021 
2022 void DofMap::dof_indices (const Elem * const elem,
2023  std::vector<dof_id_type> & di) const
2024 {
2025  // We now allow elem==NULL to request just SCALAR dofs
2026  // libmesh_assert(elem);
2027 
2028  // If we are asking for current indices on an element, it ought to
2029  // be an active element (or a Side proxy, which also thinks it's
2030  // active)
2031  libmesh_assert(!elem || elem->active());
2032 
2033  LOG_SCOPE("dof_indices()", "DofMap");
2034 
2035  // Clear the DOF indices vector
2036  di.clear();
2037 
2038  const unsigned int n_vars = this->n_variables();
2039 
2040 #ifdef DEBUG
2041  // Check that sizes match in DEBUG mode
2042  std::size_t tot_size = 0;
2043 #endif
2044 
2045  if (elem && elem->type() == TRI3SUBDIVISION)
2046  {
2047  // Subdivision surface FE require the 1-ring around elem
2048  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
2049 
2050  // Ghost subdivision elements have no real dofs
2051  if (!sd_elem->is_ghost())
2052  {
2053  // Determine the nodes contributing to element elem
2054  std::vector<const Node *> elem_nodes;
2055  MeshTools::Subdivision::find_one_ring(sd_elem, elem_nodes);
2056 
2057  // Get the dof numbers
2058  for (unsigned int v=0; v<n_vars; v++)
2059  {
2060  if(this->variable(v).type().family == SCALAR &&
2061  this->variable(v).active_on_subdomain(elem->subdomain_id()))
2062  {
2063 #ifdef DEBUG
2064  tot_size += this->variable(v).type().order;
2065 #endif
2066  std::vector<dof_id_type> di_new;
2067  this->SCALAR_dof_indices(di_new,v);
2068  di.insert( di.end(), di_new.begin(), di_new.end());
2069  }
2070  else
2071  _dof_indices(elem, elem->p_level(), di, v,
2072  &elem_nodes[0], elem_nodes.size()
2073 #ifdef DEBUG
2074  , tot_size
2075 #endif
2076  );
2077  }
2078  }
2079 
2080  return;
2081  }
2082 
2083  // Get the dof numbers
2084  for (unsigned int v=0; v<n_vars; v++)
2085  {
2086  const Variable & var = this->variable(v);
2087  if(var.type().family == SCALAR &&
2088  (!elem ||
2089  var.active_on_subdomain(elem->subdomain_id())))
2090  {
2091 #ifdef DEBUG
2092  tot_size += var.type().order;
2093 #endif
2094  std::vector<dof_id_type> di_new;
2095  this->SCALAR_dof_indices(di_new,v);
2096  di.insert( di.end(), di_new.begin(), di_new.end());
2097  }
2098  else if (elem)
2099  _dof_indices(elem, elem->p_level(), di, v, elem->get_nodes(),
2100  elem->n_nodes()
2101 #ifdef DEBUG
2102  , tot_size
2103 #endif
2104  );
2105  }
2106 
2107 #ifdef DEBUG
2108  libmesh_assert_equal_to (tot_size, di.size());
2109 #endif
2110 }
2111 
2112 
2113 void DofMap::dof_indices (const Elem * const elem,
2114  std::vector<dof_id_type> & di,
2115  const unsigned int vn,
2116  int p_level) const
2117 {
2118  // We now allow elem==NULL to request just SCALAR dofs
2119  // libmesh_assert(elem);
2120 
2121  LOG_SCOPE("dof_indices()", "DofMap");
2122 
2123  // Clear the DOF indices vector
2124  di.clear();
2125 
2126  // Use the default p refinement level?
2127  if (p_level == -12345)
2128  p_level = elem ? elem->p_level() : 0;
2129 
2130 #ifdef DEBUG
2131  // Check that sizes match in DEBUG mode
2132  std::size_t tot_size = 0;
2133 #endif
2134 
2135  if (elem && elem->type() == TRI3SUBDIVISION)
2136  {
2137  // Subdivision surface FE require the 1-ring around elem
2138  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
2139 
2140  // Ghost subdivision elements have no real dofs
2141  if (!sd_elem->is_ghost())
2142  {
2143  // Determine the nodes contributing to element elem
2144  std::vector<const Node *> elem_nodes;
2145  MeshTools::Subdivision::find_one_ring(sd_elem, elem_nodes);
2146 
2147  _dof_indices(elem, p_level, di, vn, &elem_nodes[0],
2148  elem_nodes.size()
2149 #ifdef DEBUG
2150  , tot_size
2151 #endif
2152  );
2153  }
2154 
2155  return;
2156  }
2157 
2158  const Variable & var = this->variable(vn);
2159 
2160  // Get the dof numbers
2161  if(var.type().family == SCALAR &&
2162  (!elem ||
2163  var.active_on_subdomain(elem->subdomain_id())))
2164  {
2165 #ifdef DEBUG
2166  tot_size += var.type().order;
2167 #endif
2168  std::vector<dof_id_type> di_new;
2169  this->SCALAR_dof_indices(di_new,vn);
2170  di.insert( di.end(), di_new.begin(), di_new.end());
2171  }
2172  else if (elem)
2173  _dof_indices(elem, p_level, di, vn, elem->get_nodes(),
2174  elem->n_nodes()
2175 #ifdef DEBUG
2176  , tot_size
2177 #endif
2178  );
2179 
2180 #ifdef DEBUG
2181  libmesh_assert_equal_to (tot_size, di.size());
2182 #endif
2183 }
2184 
2185 
2186 void DofMap::dof_indices (const Node * const node,
2187  std::vector<dof_id_type> & di) const
2188 {
2189  // We allow node==NULL to request just SCALAR dofs
2190  // libmesh_assert(elem);
2191 
2192  LOG_SCOPE("dof_indices(Node)", "DofMap");
2193 
2194  // Clear the DOF indices vector
2195  di.clear();
2196 
2197  const unsigned int n_vars = this->n_variables();
2198  const unsigned int sys_num = this->sys_number();
2199 
2200  // Get the dof numbers
2201  for (unsigned int v=0; v<n_vars; v++)
2202  {
2203  const Variable & var = this->variable(v);
2204  if (var.type().family == SCALAR)
2205  {
2206  std::vector<dof_id_type> di_new;
2207  this->SCALAR_dof_indices(di_new,v);
2208  di.insert( di.end(), di_new.begin(), di_new.end());
2209  }
2210  else
2211  {
2212  const int n_comp = node->n_comp(sys_num,v);
2213  for (int i=0; i != n_comp; ++i)
2214  {
2215  libmesh_assert_not_equal_to
2216  (node->dof_number(sys_num,v,i),
2218  di.push_back(node->dof_number(sys_num,v,i));
2219  }
2220  }
2221  }
2222 }
2223 
2224 
2225 void DofMap::dof_indices (const Node * const node,
2226  std::vector<dof_id_type> & di,
2227  const unsigned int vn) const
2228 {
2229  if (vn == libMesh::invalid_uint)
2230  {
2231  this->dof_indices(node, di);
2232  return;
2233  }
2234 
2235  // We allow node==NULL to request just SCALAR dofs
2236  // libmesh_assert(elem);
2237 
2238  LOG_SCOPE("dof_indices(Node)", "DofMap");
2239 
2240  // Clear the DOF indices vector
2241  di.clear();
2242 
2243  const unsigned int sys_num = this->sys_number();
2244 
2245  // Get the dof numbers
2246  const Variable & var = this->variable(vn);
2247  if (var.type().family == SCALAR)
2248  {
2249  std::vector<dof_id_type> di_new;
2250  this->SCALAR_dof_indices(di_new,vn);
2251  di.insert( di.end(), di_new.begin(), di_new.end());
2252  }
2253  else
2254  {
2255  const int n_comp = node->n_comp(sys_num,vn);
2256  for (int i=0; i != n_comp; ++i)
2257  {
2258  libmesh_assert_not_equal_to
2259  (node->dof_number(sys_num,vn,i),
2261  di.push_back(node->dof_number(sys_num,vn,i));
2262  }
2263  }
2264 }
2265 
2266 
2267 void DofMap::_dof_indices (const Elem * const elem,
2268  int p_level,
2269  std::vector<dof_id_type> & di,
2270  const unsigned int v,
2271  const Node * const * nodes,
2272  unsigned int n_nodes
2273 #ifdef DEBUG
2274  ,
2275  std::size_t & tot_size
2276 #endif
2277  ) const
2278 {
2279  // This internal function is only useful on valid elements
2280  libmesh_assert(elem);
2281 
2282  const Variable & var = this->variable(v);
2283 
2284  if (var.active_on_subdomain(elem->subdomain_id()))
2285  {
2286  const ElemType type = elem->type();
2287  const unsigned int sys_num = this->sys_number();
2288  const unsigned int dim = elem->dim();
2289 
2290  // Increase the polynomial order on p refined elements
2291  FEType fe_type = var.type();
2292  fe_type.order = static_cast<Order>(fe_type.order + p_level);
2293 
2294  const bool extra_hanging_dofs =
2296 
2297 #ifdef DEBUG
2298  // The number of dofs per element is non-static for subdivision FE
2299  if (fe_type.family == SUBDIVISION)
2300  tot_size += n_nodes;
2301  else
2302  tot_size += FEInterface::n_dofs(dim,fe_type,type);
2303 #endif
2304 
2305  // Get the node-based DOF numbers
2306  for (unsigned int n=0; n<n_nodes; n++)
2307  {
2308  const Node * node = nodes[n];
2309 
2310  // There is a potential problem with h refinement. Imagine a
2311  // quad9 that has a linear FE on it. Then, on the hanging side,
2312  // it can falsely identify a DOF at the mid-edge node. This is why
2313  // we call FEInterface instead of node->n_comp() directly.
2314  const unsigned int nc = FEInterface::n_dofs_at_node (dim,
2315  fe_type,
2316  type,
2317  n);
2318 
2319  // If this is a non-vertex on a hanging node with extra
2320  // degrees of freedom, we use the non-vertex dofs (which
2321  // come in reverse order starting from the end, to
2322  // simplify p refinement)
2323  if (extra_hanging_dofs && !elem->is_vertex(n))
2324  {
2325  const int dof_offset = node->n_comp(sys_num,v) - nc;
2326 
2327  // We should never have fewer dofs than necessary on a
2328  // node unless we're getting indices on a parent element,
2329  // and we should never need the indices on such a node
2330  if (dof_offset < 0)
2331  {
2332  libmesh_assert(!elem->active());
2333  di.resize(di.size() + nc, DofObject::invalid_id);
2334  }
2335  else
2336  for (int i=node->n_comp(sys_num,v)-1; i>=dof_offset; i--)
2337  {
2338  libmesh_assert_not_equal_to (node->dof_number(sys_num,v,i),
2340  di.push_back(node->dof_number(sys_num,v,i));
2341  }
2342  }
2343  // If this is a vertex or an element without extra hanging
2344  // dofs, our dofs come in forward order coming from the
2345  // beginning
2346  else
2347  for (unsigned int i=0; i<nc; i++)
2348  {
2349  libmesh_assert_not_equal_to (node->dof_number(sys_num,v,i),
2351  di.push_back(node->dof_number(sys_num,v,i));
2352  }
2353  }
2354 
2355  // If there are any element-based DOF numbers, get them
2356  const unsigned int nc = FEInterface::n_dofs_per_elem(dim,
2357  fe_type,
2358  type);
2359  // We should never have fewer dofs than necessary on an
2360  // element unless we're getting indices on a parent element,
2361  // and we should never need those indices
2362  if (nc != 0)
2363  {
2364  if (elem->n_systems() > sys_num &&
2365  nc <= elem->n_comp(sys_num,v))
2366  {
2367  for (unsigned int i=0; i<nc; i++)
2368  {
2369  libmesh_assert_not_equal_to (elem->dof_number(sys_num,v,i),
2371 
2372  di.push_back(elem->dof_number(sys_num,v,i));
2373  }
2374  }
2375  else
2376  {
2377  libmesh_assert(!elem->active() || fe_type.family == LAGRANGE || fe_type.family == SUBDIVISION);
2378  di.resize(di.size() + nc, DofObject::invalid_id);
2379  }
2380  }
2381  }
2382 }
2383 
2384 
2385 
2386 void DofMap::SCALAR_dof_indices (std::vector<dof_id_type> & di,
2387  const unsigned int vn,
2388 #ifdef LIBMESH_ENABLE_AMR
2389  const bool old_dofs
2390 #else
2391  const bool
2392 #endif
2393  ) const
2394 {
2395  LOG_SCOPE("SCALAR_dof_indices()", "DofMap");
2396 
2397  libmesh_assert(this->variable(vn).type().family == SCALAR);
2398 
2399 #ifdef LIBMESH_ENABLE_AMR
2400  // If we're asking for old dofs then we'd better have some
2401  if (old_dofs)
2402  libmesh_assert_greater_equal(n_old_dofs(), n_SCALAR_dofs());
2403 
2404  dof_id_type my_idx = old_dofs ?
2405  this->_first_old_scalar_df[vn] : this->_first_scalar_df[vn];
2406 #else
2407  dof_id_type my_idx = this->_first_scalar_df[vn];
2408 #endif
2409 
2410  libmesh_assert_not_equal_to(my_idx, DofObject::invalid_id);
2411 
2412  // The number of SCALAR dofs comes from the variable order
2413  const int n_dofs_vn = this->variable(vn).type().order.get_order();
2414 
2415  di.resize(n_dofs_vn);
2416  for(int i = 0; i != n_dofs_vn; ++i)
2417  di[i] = my_idx++;
2418 }
2419 
2420 
2421 
2422 bool DofMap::all_semilocal_indices (const std::vector<dof_id_type> & dof_indices_in) const
2423 {
2424  // We're all semilocal unless we find a counterexample
2425  for (std::size_t i=0; i != dof_indices_in.size(); ++i)
2426  {
2427  const dof_id_type di = dof_indices_in[i];
2428  // If it's not in the local indices
2429  if (di < this->first_dof() ||
2430  di >= this->end_dof())
2431  {
2432  // and if it's not in the ghost indices, then we're not
2433  // semilocal
2434  if (!std::binary_search(_send_list.begin(), _send_list.end(), di))
2435  return false;
2436  }
2437  }
2438 
2439  return true;
2440 }
2441 
2442 
2443 template <typename DofObjectSubclass>
2444 bool DofMap::is_evaluable(const DofObjectSubclass & obj,
2445  unsigned int var_num) const
2446 {
2447  // Everything is evaluable on a local object
2448  if (obj.processor_id() == this->processor_id())
2449  return true;
2450 
2451  std::vector<dof_id_type> di;
2452 
2453  if (var_num == libMesh::invalid_uint)
2454  this->dof_indices(&obj, di);
2455  else
2456  this->dof_indices(&obj, di, var_num);
2457 
2458  // We're all semilocal unless we find a counterexample
2459  for (std::size_t i=0; i != di.size(); ++i)
2460  {
2461  const dof_id_type dof_id = di[i];
2462  // If it's not in the local indices
2463  if (dof_id < this->first_dof() ||
2464  dof_id >= this->end_dof())
2465  {
2466  // and if it's not in the ghost indices, then we're not
2467  // evaluable
2468  if (!std::binary_search(_send_list.begin(), _send_list.end(), dof_id))
2469  return false;
2470  }
2471  }
2472 
2473  return true;
2474 }
2475 
2476 
2477 
2478 #ifdef LIBMESH_ENABLE_AMR
2479 
2480 void DofMap::old_dof_indices (const Elem * const elem,
2481  std::vector<dof_id_type> & di,
2482  const unsigned int vn) const
2483 {
2484  LOG_SCOPE("old_dof_indices()", "DofMap");
2485 
2486  libmesh_assert(elem);
2487 
2488  const ElemType type = elem->type();
2489  const unsigned int sys_num = this->sys_number();
2490  const unsigned int n_vars = this->n_variables();
2491  const unsigned int dim = elem->dim();
2492 
2493  // If we have dof indices stored on the elem, and there's no chance
2494  // that we only have those indices because we were just p refined,
2495  // then we should have old dof indices too.
2496  libmesh_assert(!elem->has_dofs(sys_num) ||
2498  elem->old_dof_object);
2499 
2500  // Clear the DOF indices vector.
2501  di.clear();
2502 
2503  // Determine the nodes contributing to element elem
2504  std::vector<const Node *> elem_nodes;
2505  if (elem->type() == TRI3SUBDIVISION)
2506  {
2507  // Subdivision surface FE require the 1-ring around elem
2508  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
2509  MeshTools::Subdivision::find_one_ring(sd_elem, elem_nodes);
2510  }
2511  else
2512  {
2513  // All other FE use only the nodes of elem itself
2514  elem_nodes.resize(elem->n_nodes(), libmesh_nullptr);
2515  for (unsigned int i=0; i<elem->n_nodes(); i++)
2516  elem_nodes[i] = elem->node_ptr(i);
2517  }
2518 
2519  // Get the dof numbers
2520  for (unsigned int v=0; v<n_vars; v++)
2521  if ((v == vn) || (vn == libMesh::invalid_uint))
2522  {
2523  if(this->variable(v).type().family == SCALAR &&
2524  (!elem ||
2525  this->variable(v).active_on_subdomain(elem->subdomain_id())))
2526  {
2527  // We asked for this variable, so add it to the vector.
2528  std::vector<dof_id_type> di_new;
2529  this->SCALAR_dof_indices(di_new,v,true);
2530  di.insert( di.end(), di_new.begin(), di_new.end());
2531  }
2532  else
2533  if (this->variable(v).active_on_subdomain(elem->subdomain_id()))
2534  { // Do this for all the variables if one was not specified
2535  // or just for the specified variable
2536 
2537  // Increase the polynomial order on p refined elements,
2538  // but make sure you get the right polynomial order for
2539  // the OLD degrees of freedom
2540  int p_adjustment = 0;
2541  if (elem->p_refinement_flag() == Elem::JUST_REFINED)
2542  {
2543  libmesh_assert_greater (elem->p_level(), 0);
2544  p_adjustment = -1;
2545  }
2546  else if (elem->p_refinement_flag() == Elem::JUST_COARSENED)
2547  {
2548  p_adjustment = 1;
2549  }
2550  FEType fe_type = this->variable_type(v);
2551  fe_type.order = static_cast<Order>(fe_type.order +
2552  elem->p_level() +
2553  p_adjustment);
2554 
2555  const bool extra_hanging_dofs =
2557 
2558  // Get the node-based DOF numbers
2559  for (std::size_t n=0; n<elem_nodes.size(); n++)
2560  {
2561  const Node * node = elem_nodes[n];
2562 
2563  // There is a potential problem with h refinement. Imagine a
2564  // quad9 that has a linear FE on it. Then, on the hanging side,
2565  // it can falsely identify a DOF at the mid-edge node. This is why
2566  // we call FEInterface instead of node->n_comp() directly.
2567  const unsigned int nc = FEInterface::n_dofs_at_node (dim,
2568  fe_type,
2569  type,
2570  n);
2572 
2573  // If this is a non-vertex on a hanging node with extra
2574  // degrees of freedom, we use the non-vertex dofs (which
2575  // come in reverse order starting from the end, to
2576  // simplify p refinement)
2577  if (extra_hanging_dofs && !elem->is_vertex(n))
2578  {
2579  const int dof_offset =
2580  node->old_dof_object->n_comp(sys_num,v) - nc;
2581 
2582  // We should never have fewer dofs than necessary on a
2583  // node unless we're getting indices on a parent element
2584  // or a just-coarsened element
2585  if (dof_offset < 0)
2586  {
2587  libmesh_assert(!elem->active() || elem->refinement_flag() ==
2589  di.resize(di.size() + nc, DofObject::invalid_id);
2590  }
2591  else
2592  for (int i=node->old_dof_object->n_comp(sys_num,v)-1;
2593  i>=dof_offset; i--)
2594  {
2595  libmesh_assert_not_equal_to (node->old_dof_object->dof_number(sys_num,v,i),
2597  di.push_back(node->old_dof_object->dof_number(sys_num,v,i));
2598  }
2599  }
2600  // If this is a vertex or an element without extra hanging
2601  // dofs, our dofs come in forward order coming from the
2602  // beginning
2603  else
2604  for (unsigned int i=0; i<nc; i++)
2605  {
2606  libmesh_assert_not_equal_to (node->old_dof_object->dof_number(sys_num,v,i),
2608  di.push_back(node->old_dof_object->dof_number(sys_num,v,i));
2609  }
2610  }
2611 
2612  // If there are any element-based DOF numbers, get them
2613  const unsigned int nc = FEInterface::n_dofs_per_elem(dim,
2614  fe_type,
2615  type);
2616 
2617  // We should never have fewer dofs than necessary on an
2618  // element unless we're getting indices on a parent element
2619  // or a just-coarsened element
2620  if (nc != 0)
2621  {
2622  if (elem->old_dof_object->n_systems() > sys_num &&
2623  nc <= elem->old_dof_object->n_comp(sys_num,v))
2624  {
2626 
2627  for (unsigned int i=0; i<nc; i++)
2628  {
2629  libmesh_assert_not_equal_to (elem->old_dof_object->dof_number(sys_num,v,i),
2631 
2632  di.push_back(elem->old_dof_object->dof_number(sys_num,v,i));
2633  }
2634  }
2635  else
2636  {
2637  libmesh_assert(!elem->active() || fe_type.family == LAGRANGE ||
2639  di.resize(di.size() + nc, DofObject::invalid_id);
2640  }
2641  }
2642  }
2643  } // end loop over variables
2644 }
2645 
2646 
2647 
2648 /*
2649  void DofMap::augment_send_list_for_projection(const MeshBase & mesh)
2650  {
2651  // Loop over the active local elements in the mesh.
2652  // If the element was just refined and its parent lives
2653  // on a different processor then we need to augment the
2654  // _send_list with the parent's DOF indices so that we
2655  // can access the parent's data for computing solution
2656  // projections, etc...
2657 
2658  // The DOF indices for the parent
2659  std::vector<dof_id_type> di;
2660 
2661  // Flag telling us if we need to re-sort the send_list.
2662  // Note we won't need to re-sort unless a child with a
2663  // parent belonging to a different processor is found
2664  bool needs_sorting = false;
2665 
2666 
2667  MeshBase::const_element_iterator elem_it = mesh.active_local_elements_begin();
2668  const MeshBase::const_element_iterator elem_end = mesh.active_local_elements_end();
2669 
2670  for ( ; elem_it != elem_end; ++elem_it)
2671  {
2672  const Elem * const elem = *elem_it;
2673 
2674  // We only need to consider the children that
2675  // were just refined
2676  if (elem->refinement_flag() != Elem::JUST_REFINED) continue;
2677 
2678  const Elem * const parent = elem->parent();
2679 
2680  // If the parent lives on another processor
2681  // than the child
2682  if (parent != libmesh_nullptr)
2683  if (parent->processor_id() != elem->processor_id())
2684  {
2685  // Get the DOF indices for the parent
2686  this->dof_indices (parent, di);
2687 
2688  // Insert the DOF indices into the send list
2689  _send_list.insert (_send_list.end(),
2690  di.begin(), di.end());
2691 
2692  // We will need to re-sort the send list
2693  needs_sorting = true;
2694  }
2695  }
2696 
2697  // The send-list might need to be sorted again.
2698  if (needs_sorting) this->sort_send_list ();
2699  }
2700 */
2701 
2702 #endif // LIBMESH_ENABLE_AMR
2703 
2704 
2705 #ifdef LIBMESH_ENABLE_CONSTRAINTS
2706 
2707 void DofMap::find_connected_dofs (std::vector<dof_id_type> & elem_dofs) const
2708 {
2709  typedef std::set<dof_id_type> RCSet;
2710 
2711  // First insert the DOFS we already depend on into the set.
2712  RCSet dof_set (elem_dofs.begin(), elem_dofs.end());
2713 
2714  bool done = true;
2715 
2716  // Next insert any dofs those might be constrained in terms
2717  // of. Note that in this case we may not be done: Those may
2718  // in turn depend on others. So, we need to repeat this process
2719  // in that case until the system depends only on unconstrained
2720  // degrees of freedom.
2721  for (std::size_t i=0; i<elem_dofs.size(); i++)
2722  if (this->is_constrained_dof(elem_dofs[i]))
2723  {
2724  // If the DOF is constrained
2725  DofConstraints::const_iterator
2726  pos = _dof_constraints.find(elem_dofs[i]);
2727 
2728  libmesh_assert (pos != _dof_constraints.end());
2729 
2730  const DofConstraintRow & constraint_row = pos->second;
2731 
2732  // adaptive p refinement currently gives us lots of empty constraint
2733  // rows - we should optimize those DoFs away in the future. [RHS]
2734  //libmesh_assert (!constraint_row.empty());
2735 
2736  DofConstraintRow::const_iterator it = constraint_row.begin();
2737  DofConstraintRow::const_iterator it_end = constraint_row.end();
2738 
2739 
2740  // Add the DOFs this dof is constrained in terms of.
2741  // note that these dofs might also be constrained, so
2742  // we will need to call this function recursively.
2743  for ( ; it != it_end; ++it)
2744  if (!dof_set.count (it->first))
2745  {
2746  dof_set.insert (it->first);
2747  done = false;
2748  }
2749  }
2750 
2751 
2752  // If not done then we need to do more work
2753  // (obviously :-) )!
2754  if (!done)
2755  {
2756  // Fill the vector with the contents of the set
2757  elem_dofs.clear();
2758  elem_dofs.insert (elem_dofs.end(),
2759  dof_set.begin(), dof_set.end());
2760 
2761 
2762  // May need to do this recursively. It is possible
2763  // that we just replaced a constrained DOF with another
2764  // constrained DOF.
2765  this->find_connected_dofs (elem_dofs);
2766 
2767  } // end if (!done)
2768 }
2769 
2770 #endif // LIBMESH_ENABLE_CONSTRAINTS
2771 
2772 
2773 
2774 void DofMap::print_info(std::ostream & os) const
2775 {
2776  os << this->get_info();
2777 }
2778 
2779 
2780 
2781 std::string DofMap::get_info() const
2782 {
2783  std::ostringstream os;
2784 
2785  // If we didn't calculate the exact sparsity pattern, the threaded
2786  // sparsity pattern assembly may have just given us an upper bound
2787  // on sparsity.
2788  const char * may_equal = " <= ";
2789 
2790  // If we calculated the exact sparsity pattern, then we can report
2791  // exact bandwidth figures:
2792  std::vector<SparseMatrix<Number> * >::const_iterator
2793  pos = _matrices.begin(),
2794  end = _matrices.end();
2795 
2796  for (; pos != end; ++pos)
2797  if ((*pos)->need_full_sparsity_pattern())
2798  may_equal = " = ";
2799 
2800  dof_id_type max_n_nz = 0, max_n_oz = 0;
2801  long double avg_n_nz = 0, avg_n_oz = 0;
2802 
2803  if (_n_nz)
2804  {
2805  for (std::size_t i = 0; i != _n_nz->size(); ++i)
2806  {
2807  max_n_nz = std::max(max_n_nz, (*_n_nz)[i]);
2808  avg_n_nz += (*_n_nz)[i];
2809  }
2810 
2811  std::size_t n_nz_size = _n_nz->size();
2812 
2813  this->comm().max(max_n_nz);
2814  this->comm().sum(avg_n_nz);
2815  this->comm().sum(n_nz_size);
2816 
2817  avg_n_nz /= std::max(n_nz_size,std::size_t(1));
2818 
2820 
2821  for (std::size_t i = 0; i != (*_n_oz).size(); ++i)
2822  {
2823  max_n_oz = std::max(max_n_oz, (*_n_oz)[i]);
2824  avg_n_oz += (*_n_oz)[i];
2825  }
2826 
2827  std::size_t n_oz_size = _n_oz->size();
2828 
2829  this->comm().max(max_n_oz);
2830  this->comm().sum(avg_n_oz);
2831  this->comm().sum(n_oz_size);
2832 
2833  avg_n_oz /= std::max(n_oz_size,std::size_t(1));
2834  }
2835 
2836  os << " DofMap Sparsity\n Average On-Processor Bandwidth"
2837  << may_equal << avg_n_nz << '\n';
2838 
2839  os << " Average Off-Processor Bandwidth"
2840  << may_equal << avg_n_oz << '\n';
2841 
2842  os << " Maximum On-Processor Bandwidth"
2843  << may_equal << max_n_nz << '\n';
2844 
2845  os << " Maximum Off-Processor Bandwidth"
2846  << may_equal << max_n_oz << std::endl;
2847 
2848 #ifdef LIBMESH_ENABLE_CONSTRAINTS
2849 
2850  std::size_t n_constraints = 0, max_constraint_length = 0,
2851  n_rhss = 0;
2852  long double avg_constraint_length = 0.;
2853 
2854  for (DofConstraints::const_iterator it=_dof_constraints.begin();
2855  it != _dof_constraints.end(); ++it)
2856  {
2857  // Only count local constraints, then sum later
2858  const dof_id_type constrained_dof = it->first;
2859  if (constrained_dof < this->first_dof() ||
2860  constrained_dof >= this->end_dof())
2861  continue;
2862 
2863  const DofConstraintRow & row = it->second;
2864  std::size_t rowsize = row.size();
2865 
2866  max_constraint_length = std::max(max_constraint_length,
2867  rowsize);
2868  avg_constraint_length += rowsize;
2869  n_constraints++;
2870 
2871  if (_primal_constraint_values.count(constrained_dof))
2872  n_rhss++;
2873  }
2874 
2875  this->comm().sum(n_constraints);
2876  this->comm().sum(n_rhss);
2877  this->comm().sum(avg_constraint_length);
2878  this->comm().max(max_constraint_length);
2879 
2880  os << " DofMap Constraints\n Number of DoF Constraints = "
2881  << n_constraints;
2882  if (n_rhss)
2883  os << '\n'
2884  << " Number of Heterogenous Constraints= " << n_rhss;
2885  if (n_constraints)
2886  {
2887  avg_constraint_length /= n_constraints;
2888 
2889  os << '\n'
2890  << " Average DoF Constraint Length= " << avg_constraint_length;
2891  }
2892 
2893 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2894  std::size_t n_node_constraints = 0, max_node_constraint_length = 0,
2895  n_node_rhss = 0;
2896  long double avg_node_constraint_length = 0.;
2897 
2898  for (NodeConstraints::const_iterator it=_node_constraints.begin();
2899  it != _node_constraints.end(); ++it)
2900  {
2901  // Only count local constraints, then sum later
2902  const Node * node = it->first;
2903  if (node->processor_id() != this->processor_id())
2904  continue;
2905 
2906  const NodeConstraintRow & row = it->second.first;
2907  std::size_t rowsize = row.size();
2908 
2909  max_node_constraint_length = std::max(max_node_constraint_length,
2910  rowsize);
2911  avg_node_constraint_length += rowsize;
2912  n_node_constraints++;
2913 
2914  if (it->second.second != Point(0))
2915  n_node_rhss++;
2916  }
2917 
2918  this->comm().sum(n_node_constraints);
2919  this->comm().sum(n_node_rhss);
2920  this->comm().sum(avg_node_constraint_length);
2921  this->comm().max(max_node_constraint_length);
2922 
2923  os << "\n Number of Node Constraints = " << n_node_constraints;
2924  if (n_node_rhss)
2925  os << '\n'
2926  << " Number of Heterogenous Node Constraints= " << n_node_rhss;
2927  if (n_node_constraints)
2928  {
2929  avg_node_constraint_length /= n_node_constraints;
2930  os << "\n Maximum Node Constraint Length= " << max_node_constraint_length
2931  << '\n'
2932  << " Average Node Constraint Length= " << avg_node_constraint_length;
2933  }
2934 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2935 
2936  os << std::endl;
2937 
2938 #endif // LIBMESH_ENABLE_CONSTRAINTS
2939 
2940  return os.str();
2941 }
2942 
2943 
2944 template bool DofMap::is_evaluable<Elem>(const Elem &, unsigned int) const;
2945 template bool DofMap::is_evaluable<Node>(const Node &, unsigned int) const;
2946 
2947 } // namespace libMesh
std::vector< VariableGroup > _variable_groups
Definition: dof_map.h:1384
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:2274
FEFamily family
Definition: fe_type.h:206
int get_order() const
Definition: fe_type.h:77
bool _implicit_neighbor_dofs_initialized
Definition: dof_map.h:1612
void distribute_local_dofs_node_major(dof_id_type &next_free_dof, MeshBase &mesh)
Definition: dof_map.C:1199
const FEType & type() const
Definition: variable.h:119
const unsigned int _sys_number
Definition: dof_map.h:1389
bool _implicit_neighbor_dofs
Definition: dof_map.h:1613
std::string get_info() const
Definition: dof_map.C:2781
void print_info(std::ostream &os=libMesh::out) const
Definition: dof_map.C:2774
UniquePtr< DefaultCoupling > _default_coupling
Definition: dof_map.h:1463
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:1984
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:1440
unsigned int p_level() const
Definition: elem.h:2149
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1668
bool is_evaluable(const DofObjectSubclass &obj, unsigned int var_num=libMesh::invalid_uint) const
Definition: dof_map.C:2444
void set_implicit_neighbor_dofs(bool implicit_neighbor_dofs)
Definition: dof_map.C:1775
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:1423
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:303
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:2267
const class libmesh_nullptr_t libmesh_nullptr
AugmentSparsityPattern * _augment_sparsity_pattern
Definition: dof_map.h:1428
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:1481
std::vector< dof_id_type > _end_old_df
Definition: dof_map.h:1548
OrderWrapper order
Definition: fe_type.h:200
void add_coupling_functor(GhostingFunctor &coupling_functor)
Definition: dof_map.C:1898
IterBase * end
void local_variable_indices(std::vector< dof_id_type > &idx, const MeshBase &mesh, unsigned int var_num) const
Definition: dof_map.C:1110
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:821
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:2480
const_iterator end() const
const Variable & variable(const unsigned int c) const
Definition: dof_map.h:1638
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:1471
dof_id_type n_dofs() const
Definition: dof_map.h:506
virtual void zero()=0
const VariableGroup & variable_group(const unsigned int c) const
Definition: dof_map.h:1628
long double max(long double a, double b)
std::vector< dof_id_type > _first_old_df
Definition: dof_map.h:1543
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:1417
dof_id_type n_dofs_on_processor(const processor_id_type proc) const
Definition: dof_map.h:522
AdjointDofConstraintValues _adjoint_constraint_values
Definition: dof_map.h:1567
libmesh_assert(j)
std::vector< dof_id_type > * _n_nz
Definition: dof_map.h:1514
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
AugmentSendList * _augment_send_list
Definition: dof_map.h:1445
std::vector< dof_id_type > * _n_oz
Definition: dof_map.h:1520
bool need_full_sparsity_pattern
Definition: dof_map.h:1500
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:1917
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:1506
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:1705
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:1574
std::vector< dof_id_type > _first_old_scalar_df
Definition: dof_map.h:1554
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:1334
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:1433
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1658
dof_id_type _n_SCALAR_dofs
Definition: dof_map.h:1531
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:1926
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:346
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:215
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:134
virtual bool need_full_sparsity_pattern() const
void clear()
Definition: dof_map.C:842
const Node & node_ref(const unsigned int i) const
Definition: elem.h:1680
A surface shell element used in mechanics calculations.
static const dof_id_type invalid_id
Definition: dof_object.h:335
void distribute_dofs(MeshBase &)
Definition: dof_map.C:922
RefinementState p_refinement_flag() const
Definition: elem.h:2236
DofConstraints _dof_constraints
Definition: dof_map.h:1563
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:1210
unsigned int n_variable_groups() const
Definition: dof_map.h:465
OStreamProxy err(std::cerr)
subdomain_id_type subdomain_id() const
Definition: elem.h:1733
std::vector< DirichletBoundaries * > _adjoint_dirichlet_boundaries
Definition: dof_map.h:1603
unsigned int sys_number() const
Definition: dof_map.h:1620
bool is_prepared() const
Definition: mesh_base.h:127
void * _extra_send_list_context
Definition: dof_map.h:1455
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:1499
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:1401
dof_id_type end_dof() const
Definition: dof_map.h:572
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:473
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:1597
void SCALAR_dof_indices(std::vector< dof_id_type > &di, const unsigned int vn, const bool old_dofs=false) const
Definition: dof_map.C:2386
void remove_coupling_functor(GhostingFunctor &coupling_functor)
Definition: dof_map.C:1907
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:1163
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:1565
std::vector< Variable > _variables
Definition: dof_map.h:1379
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
DofConstraints _stashed_dof_constraints
Definition: dof_map.h:1563
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:1538
void clear_sparsity()
Definition: dof_map.C:1876
bool all_semilocal_indices(const std::vector< dof_id_type > &dof_indices) const
Definition: dof_map.C:2422
RefinementState refinement_flag() const
Definition: elem.h:2220
void parallel_reduce(const Range &range, Body &body)
Definition: threads_none.h:101
std::set< GhostingFunctor * > _coupling_functors
Definition: dof_map.h:1494
virtual numeric_index_type first_local_index() const =0
UniquePtr< PeriodicBoundaries > _periodic_boundaries
Definition: dof_map.h:1583
std::vector< dof_id_type > _end_df
Definition: dof_map.h:1411
dof_id_type n_SCALAR_dofs() const
Definition: dof_map.h:511
std::vector< dof_id_type > _first_df
Definition: dof_map.h:1406
void(* _extra_send_list_function)(std::vector< dof_id_type > &, void *)
Definition: dof_map.h:1450
DofObject *(DofMap::* dofobject_accessor)(MeshBase &mesh, dof_id_type i) const
Definition: dof_map.h:1262
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:1394
void find_connected_dofs(std::vector< dof_id_type > &elem_dofs) const
Definition: dof_map.C:2707
void prepare_send_list()
Definition: dof_map.C:1739
Real size() const
Definition: type_vector.h:877
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:1839
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:1782
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:1936
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:1650
void distribute_local_dofs_var_major(dof_id_type &next_free_dof, MeshBase &mesh)
Definition: dof_map.C:1353
dof_id_type first_dof() const
Definition: dof_map.h:534
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:743
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:1525
processor_id_type processor_id() const
Definition: dof_object.h:686
NodeConstraints _node_constraints
Definition: dof_map.h:1574
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:2022