libMesh
dof_map.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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
21 #include "libmesh/dof_map.h"
22 
23 // libMesh includes
24 #include "libmesh/coupling_matrix.h"
25 #include "libmesh/default_coupling.h"
26 #include "libmesh/dense_matrix.h"
27 #include "libmesh/dense_vector_base.h"
28 #include "libmesh/dirichlet_boundaries.h"
29 #include "libmesh/enum_to_string.h"
30 #include "libmesh/fe_type.h"
31 #include "libmesh/fe_base.h" // FEBase::build() for continuity test
32 #include "libmesh/ghosting_functor.h"
33 #include "libmesh/int_range.h"
34 #include "libmesh/mesh_base.h"
35 #include "libmesh/mesh_tools.h"
36 #include "libmesh/numeric_vector.h"
37 #include "libmesh/periodic_boundary_base.h"
38 #include "libmesh/periodic_boundaries.h"
39 #include "libmesh/sparse_matrix.h"
40 #include "libmesh/sparsity_pattern.h"
41 #include "libmesh/threads.h"
42 #include "libmesh/static_condensation_dof_map.h"
43 
44 // TIMPI includes
46 #include "timpi/parallel_sync.h"
47 
48 // C++ Includes
49 #include <algorithm> // for std::fill, std::equal_range, std::max, std::lower_bound, etc.
50 #include <memory>
51 #include <set>
52 #include <sstream>
53 #include <unordered_map>
54 
55 namespace libMesh
56 {
57 
58 // ------------------------------------------------------------
59 // DofMap member functions
60 std::unique_ptr<SparsityPattern::Build>
62  const bool calculate_constrained,
63  const bool use_condensed_system) const
64 {
65  libmesh_assert (mesh.is_prepared());
66 
67  LOG_SCOPE("build_sparsity()", "DofMap");
68 
69  // Compute the sparsity structure of the global matrix. This can be
70  // fed into a PetscMatrixBase to allocate exactly the number of nonzeros
71  // necessary to store the matrix. This algorithm should be linear
72  // in the (# of elements)*(# nodes per element)
73 
74  // We can be more efficient in the threaded sparsity pattern assembly
75  // if we don't need the exact pattern. For some sparse matrix formats
76  // a good upper bound will suffice.
77 
78  // See if we need to include sparsity pattern entries for coupling
79  // between neighbor dofs
80  bool implicit_neighbor_dofs = this->use_coupled_neighbor_dofs(mesh);
81 
82  const StaticCondensationDofMap * sc = nullptr;
83  if (use_condensed_system)
84  {
86  sc = _sc.get();
87  }
88 
89  // We can compute the sparsity pattern in parallel on multiple
90  // threads. The goal is for each thread to compute the full sparsity
91  // pattern for a subset of elements. These sparsity patterns can
92  // be efficiently merged in the SparsityPattern::Build::join()
93  // method, especially if there is not too much overlap between them.
94  // Even better, if the full sparsity pattern is not needed then
95  // the number of nonzeros per row can be estimated from the
96  // sparsity patterns created on each thread.
97  auto sp = std::make_unique<SparsityPattern::Build>
98  (*this,
99  this->_dof_coupling,
100  this->_coupling_functors,
101  implicit_neighbor_dofs,
103  calculate_constrained,
104  sc);
105 
106  Threads::parallel_reduce (ConstElemRange (mesh.active_local_elements_begin(),
107  mesh.active_local_elements_end()), *sp);
108 
109  sp->parallel_sync();
110 
111  libmesh_assert_equal_to (sp->get_sparsity_pattern().size(), this->n_local_dofs());
112 
113  // Check to see if we have any extra stuff to add to the sparsity_pattern
115  {
117  {
118  libmesh_here();
119  libMesh::out << "WARNING: You have specified both an extra sparsity function and object.\n"
120  << " Are you sure this is what you meant to do??"
121  << std::endl;
122  }
123 
124  sp->apply_extra_sparsity_function(_extra_sparsity_function,
126  }
127 
129  sp->apply_extra_sparsity_object(*_augment_sparsity_pattern);
130 
131  return sp;
132 }
133 
134 
135 
136 DofMap::DofMap(const unsigned int number,
137  MeshBase & mesh) :
138  DofMapBase (mesh.comm()),
139  _dof_coupling(nullptr),
140  _error_on_constraint_loop(false),
141  _constrained_sparsity_construction(false),
142  _variables(),
143  _variable_groups(),
144  _variable_group_numbers(),
145  _sys_number(number),
146  _mesh(mesh),
147  _matrices(),
148  _first_scalar_df(),
149  _send_list(),
150  _augment_sparsity_pattern(nullptr),
151  _extra_sparsity_function(nullptr),
152  _extra_sparsity_context(nullptr),
153  _augment_send_list(nullptr),
154  _extra_send_list_function(nullptr),
155  _extra_send_list_context(nullptr),
156  _default_coupling(std::make_unique<DefaultCoupling>()),
157  _default_evaluating(std::make_unique<DefaultCoupling>()),
158  need_full_sparsity_pattern(false),
159  _n_SCALAR_dofs(0)
160 #ifdef LIBMESH_ENABLE_AMR
161  , _first_old_scalar_df()
162 #endif
163 #ifdef LIBMESH_ENABLE_CONSTRAINTS
164  , _dof_constraints()
165  , _stashed_dof_constraints()
166  , _primal_constraint_values()
167  , _adjoint_constraint_values()
168 #endif
169 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
170  , _node_constraints()
171 #endif
172 #ifdef LIBMESH_ENABLE_PERIODIC
173  , _periodic_boundaries(std::make_unique<PeriodicBoundaries>())
174 #endif
175 #ifdef LIBMESH_ENABLE_DIRICHLET
176  , _dirichlet_boundaries(std::make_unique<DirichletBoundaries>())
177  , _adjoint_dirichlet_boundaries()
178 #endif
179  , _implicit_neighbor_dofs_initialized(false),
180  _implicit_neighbor_dofs(false),
181  _verify_dirichlet_bc_consistency(true),
182  _sc(nullptr)
183 {
184  _matrices.clear();
185 
186  _default_coupling->set_mesh(&_mesh);
187  _default_evaluating->set_mesh(&_mesh);
188  _default_evaluating->set_n_levels(1);
189 
190 #ifdef LIBMESH_ENABLE_PERIODIC
191  _default_coupling->set_periodic_boundaries(_periodic_boundaries.get());
192  _default_evaluating->set_periodic_boundaries(_periodic_boundaries.get());
193 #endif
194 
197 }
198 
199 
200 
201 // Destructor
203 {
204  this->clear();
205 
206  // clear() resets all but the default DofMap-based functors. We
207  // need to remove those from the mesh too before we die.
210 }
211 
212 
213 #ifdef LIBMESH_ENABLE_PERIODIC
214 
215 bool DofMap::is_periodic_boundary (const boundary_id_type boundaryid) const
216 {
217  if (_periodic_boundaries->count(boundaryid) != 0)
218  return true;
219 
220  return false;
221 }
222 
223 #endif
224 
225 
226 
227 // void DofMap::add_variable (const Variable & var)
228 // {
229 // libmesh_not_implemented();
230 // _variables.push_back (var);
231 // }
232 
233 
234 void DofMap::set_error_on_cyclic_constraint(bool error_on_cyclic_constraint)
235 {
236  // This function will eventually be officially libmesh_deprecated();
237  // Call DofMap::set_error_on_constraint_loop() instead.
238  set_error_on_constraint_loop(error_on_cyclic_constraint);
239 }
240 
241 void DofMap::set_error_on_constraint_loop(bool error_on_constraint_loop)
242 {
243  _error_on_constraint_loop = error_on_constraint_loop;
244 }
245 
246 
247 
249 {
250  // Ensure that we are not duplicating an existing entry in _variable_groups
251  if (std::find(_variable_groups.begin(), _variable_groups.end(), var_group) == _variable_groups.end())
252  {
253  const unsigned int vg = cast_int<unsigned int>(_variable_groups.size());
254 
255  _variable_groups.push_back(std::move(var_group));
256 
257  VariableGroup & new_var_group = _variable_groups.back();
258 
259  for (auto var : make_range(new_var_group.n_variables()))
260  {
261  auto var_instance = new_var_group(var);
262  const auto vn = var_instance.number();
263  _variables.push_back (std::move(var_instance));
264  _variable_group_numbers.push_back (vg);
265  _var_to_vg.emplace(vn, vg);
266  }
267  }
268  // End if check for var_group in _variable_groups
269 
270 }
271 
272 
273 
275 {
276  parallel_object_only();
277 
278  // We shouldn't be trying to re-attach the same matrices repeatedly
279  libmesh_assert (std::find(_matrices.begin(), _matrices.end(),
280  &matrix) == _matrices.end());
281 
282  _matrices.push_back(&matrix);
283 
284  this->update_sparsity_pattern(matrix);
285 
286  if (matrix.need_full_sparsity_pattern())
288 }
289 
290 
291 
293 {
295  (!_sp->get_n_nz().empty() ||
296  !_sp->get_n_oz().empty());
297  this->comm().max(computed_sparsity_already);
299 }
300 
301 
302 
304 {
305  matrix.attach_dof_map (*this);
306 
307  // If we've already computed sparsity, then it's too late
308  // to wait for "compute_sparsity" to help with sparse matrix
309  // initialization, and we need to handle this matrix individually
310  if (this->computed_sparsity_already())
311  {
312  libmesh_assert(_sp.get());
313 
314  if (matrix.need_full_sparsity_pattern())
315  {
316  // We'd better have already computed the full sparsity
317  // pattern if we need it here
319 
320  matrix.update_sparsity_pattern (_sp->get_sparsity_pattern());
321  }
322 
323  matrix.attach_sparsity_pattern(*_sp);
324  }
325 }
326 
327 
328 
330 {
331  return (std::find(_matrices.begin(), _matrices.end(),
332  &matrix) != _matrices.end());
333 }
334 
335 
336 
338 {
339  return mesh.node_ptr(i);
340 }
341 
342 
343 
345 {
346  return mesh.elem_ptr(i);
347 }
348 
349 
350 
351 template <typename iterator_type>
352 void DofMap::set_nonlocal_dof_objects(iterator_type objects_begin,
353  iterator_type objects_end,
354  MeshBase & mesh,
355  dofobject_accessor objects)
356 {
357  // This function must be run on all processors at once
358  parallel_object_only();
359 
360  // First, iterate over local objects to find out how many
361  // are on each processor
362  std::unordered_map<processor_id_type, dof_id_type> ghost_objects_from_proc;
363 
364  iterator_type it = objects_begin;
365 
366  for (; it != objects_end; ++it)
367  {
368  DofObject * obj = *it;
369 
370  if (obj)
371  {
372  processor_id_type obj_procid = obj->processor_id();
373  // We'd better be completely partitioned by now
374  libmesh_assert_not_equal_to (obj_procid, DofObject::invalid_processor_id);
375  ghost_objects_from_proc[obj_procid]++;
376  }
377  }
378 
379  // Request sets to send to each processor
380  std::map<processor_id_type, std::vector<dof_id_type>>
381  requested_ids;
382 
383  // We know how many of our objects live on each processor, so
384  // reserve() space for requests from each.
385  for (auto [p, size] : ghost_objects_from_proc)
386  {
387  if (p != this->processor_id())
388  requested_ids[p].reserve(size);
389  }
390 
391  for (it = objects_begin; it != objects_end; ++it)
392  {
393  DofObject * obj = *it;
395  requested_ids[obj->processor_id()].push_back(obj->id());
396  }
397 #ifdef DEBUG
398  for (auto p : make_range(this->n_processors()))
399  {
400  if (ghost_objects_from_proc.count(p))
401  libmesh_assert_equal_to (requested_ids[p].size(), ghost_objects_from_proc[p]);
402  else
403  libmesh_assert(!requested_ids.count(p));
404  }
405 #endif
406 
407  typedef std::vector<dof_id_type> datum;
408 
409  auto gather_functor =
410  [this, &mesh, &objects]
412  const std::vector<dof_id_type> & ids,
413  std::vector<datum> & data)
414  {
415  // Fill those requests
416  const unsigned int
417  sys_num = this->sys_number(),
418  n_var_groups = this->n_variable_groups();
419 
420  const std::size_t query_size = ids.size();
421 
422  data.resize(query_size);
423  for (auto & d : data)
424  d.resize(2 * n_var_groups);
425 
426  for (std::size_t i=0; i != query_size; ++i)
427  {
428  DofObject * requested = (this->*objects)(mesh, ids[i]);
429  libmesh_assert(requested);
430  libmesh_assert_equal_to (requested->processor_id(), this->processor_id());
431  libmesh_assert_equal_to (requested->n_var_groups(sys_num), n_var_groups);
432  for (unsigned int vg=0; vg != n_var_groups; ++vg)
433  {
434  unsigned int n_comp_g =
435  requested->n_comp_group(sys_num, vg);
436  data[i][vg] = n_comp_g;
437  dof_id_type my_first_dof = n_comp_g ?
438  requested->vg_dof_base(sys_num, vg) : 0;
439  libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
440  data[i][n_var_groups+vg] = my_first_dof;
441  }
442  }
443  };
444 
445  auto action_functor =
446  [this, &mesh, &objects]
447  (processor_id_type libmesh_dbg_var(pid),
448  const std::vector<dof_id_type> & ids,
449  const std::vector<datum> & data)
450  {
451  const unsigned int
452  sys_num = this->sys_number(),
453  n_var_groups = this->n_variable_groups();
454 
455  // Copy the id changes we've now been informed of
456  for (auto i : index_range(ids))
457  {
458  DofObject * requested = (this->*objects)(mesh, ids[i]);
459  libmesh_assert(requested);
460  libmesh_assert_equal_to (requested->processor_id(), pid);
461  for (unsigned int vg=0; vg != n_var_groups; ++vg)
462  {
463  unsigned int n_comp_g =
464  cast_int<unsigned int>(data[i][vg]);
465  requested->set_n_comp_group(sys_num, vg, n_comp_g);
466  if (n_comp_g)
467  {
468  dof_id_type my_first_dof = data[i][n_var_groups+vg];
469  libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
470  requested->set_vg_dof_base
471  (sys_num, vg, my_first_dof);
472  }
473  }
474  }
475  };
476 
477  datum * ex = nullptr;
478  Parallel::pull_parallel_vector_data
479  (this->comm(), requested_ids, gather_functor, action_functor, ex);
480 
481 #ifdef DEBUG
482  // Double check for invalid dofs
483  for (it = objects_begin; it != objects_end; ++it)
484  {
485  DofObject * obj = *it;
486  libmesh_assert (obj);
487  unsigned int num_variables = obj->n_vars(this->sys_number());
488  for (unsigned int v=0; v != num_variables; ++v)
489  {
490  unsigned int n_comp =
491  obj->n_comp(this->sys_number(), v);
492  dof_id_type my_first_dof = n_comp ?
493  obj->dof_number(this->sys_number(), v, 0) : 0;
494  libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
495  }
496  }
497 #endif
498 }
499 
500 
501 
502 void DofMap::reinit
504  const std::map<const Node *, std::set<subdomain_id_type>> &
505  constraining_subdomains)
506 {
508 
509  LOG_SCOPE("reinit()", "DofMap");
510 
511  // This is the common case and we want to optimize for it
512  const bool constraining_subdomains_empty =
513  constraining_subdomains.empty();
514 
515  // We ought to reconfigure our default coupling functor.
516  //
517  // The user might have removed it from our coupling functors set,
518  // but if so, who cares, this reconfiguration is cheap.
519 
520  // Avoid calling set_dof_coupling() with an empty/non-nullptr
521  // _dof_coupling matrix which may happen when there are actually no
522  // variables on the system.
523  if (this->_dof_coupling && this->_dof_coupling->empty() && !this->n_variables())
524  this->_dof_coupling = nullptr;
525  _default_coupling->set_dof_coupling(this->_dof_coupling);
526 
527  // By default we may want 0 or 1 levels of coupling
528  unsigned int standard_n_levels =
529  this->use_coupled_neighbor_dofs(mesh);
530  _default_coupling->set_n_levels
531  (std::max(_default_coupling->n_levels(), standard_n_levels));
532 
533  // But we *don't* want to restrict to a CouplingMatrix unless the
534  // user does so manually; the original libMesh behavior was to put
535  // ghost indices on the send_list regardless of variable.
536  //_default_evaluating->set_dof_coupling(this->_dof_coupling);
537 
538  const unsigned int
539  sys_num = this->sys_number(),
540  n_var_groups = this->n_variable_groups();
541 
542  // The DofObjects need to know how many variable groups we have, and
543  // how many variables there are in each group.
544  std::vector<unsigned int> n_vars_per_group; n_vars_per_group.reserve (n_var_groups);
545 
546  for (unsigned int vg=0; vg<n_var_groups; vg++)
547  n_vars_per_group.push_back (this->variable_group(vg).n_variables());
548 
549 #ifdef LIBMESH_ENABLE_AMR
550 
551  //------------------------------------------------------------
552  // Clear the old_dof_objects for all the nodes
553  // and elements so that we can overwrite them
554  for (auto & node : mesh.node_ptr_range())
555  {
556  node->clear_old_dof_object();
557  libmesh_assert (!node->get_old_dof_object());
558  }
559 
560  for (auto & elem : mesh.element_ptr_range())
561  {
562  elem->clear_old_dof_object();
563  libmesh_assert (!elem->get_old_dof_object());
564  }
565 
566 
567  //------------------------------------------------------------
568  // Set the old_dof_objects for the elements that
569  // weren't just created, if these old dof objects
570  // had variables
571  for (auto & elem : mesh.element_ptr_range())
572  {
573  // Skip the elements that were just refined
574  if (elem->refinement_flag() == Elem::JUST_REFINED)
575  continue;
576 
577  for (Node & node : elem->node_ref_range())
578  if (node.get_old_dof_object() == nullptr)
579  if (node.has_dofs(sys_num))
580  node.set_old_dof_object();
581 
582  libmesh_assert (!elem->get_old_dof_object());
583 
584  if (elem->has_dofs(sys_num))
585  elem->set_old_dof_object();
586  }
587 
588 #endif // #ifdef LIBMESH_ENABLE_AMR
589 
590 
591  //------------------------------------------------------------
592  // Then set the number of variables for each \p DofObject
593  // equal to n_variables() for this system. This will
594  // handle new \p DofObjects that may have just been created
595 
596  // All the nodes
597  for (auto & node : mesh.node_ptr_range())
598  node->set_n_vars_per_group(sys_num, n_vars_per_group);
599 
600  // All the elements
601  for (auto & elem : mesh.element_ptr_range())
602  elem->set_n_vars_per_group(sys_num, n_vars_per_group);
603 
604  // Zero _n_SCALAR_dofs, it will be updated below.
605  this->_n_SCALAR_dofs = 0;
606 
607  //------------------------------------------------------------
608  // Next allocate space for the DOF indices
609  for (unsigned int vg=0; vg<n_var_groups; vg++)
610  {
611  const VariableGroup & vg_description = this->variable_group(vg);
612 
613  const unsigned int n_var_in_group = vg_description.n_variables();
614  const FEType & base_fe_type = vg_description.type();
615 
616  const bool add_p_level =
617 #ifdef LIBMESH_ENABLE_AMR
618  !_dont_p_refine.count(vg);
619 #else
620  false;
621 #endif
622 
623  // Don't need to loop over elements for a SCALAR variable
624  // Just increment _n_SCALAR_dofs
625  if (base_fe_type.family == SCALAR)
626  {
627  this->_n_SCALAR_dofs += base_fe_type.order.get_order()*n_var_in_group;
628  continue;
629  }
630 
631  // This should be constant even on p-refined elements
632  const bool extra_hanging_dofs =
633  FEInterface::extra_hanging_dofs(base_fe_type);
634 
635  // For all the active elements, count vertex degrees of freedom.
636  for (auto & elem : mesh.active_element_ptr_range())
637  {
638  libmesh_assert(elem);
639 
640  // Only number dofs connected to active elements on this
641  // processor and only for variables which are active on on
642  // this element's subdomain or which are active on the
643  // subdomain of a node constrained by this node.
644  const bool active_on_elem =
645  vg_description.active_on_subdomain(elem->subdomain_id());
646 
647  // If there's no way we're active on this element then we're
648  // done
649  if (!active_on_elem && constraining_subdomains_empty)
650  continue;
651 
652  FEType fe_type = base_fe_type;
653 
654  const ElemType type = elem->type();
655 
656  libmesh_error_msg_if(base_fe_type.order.get_order() >
657  int(FEInterface::max_order(base_fe_type,type)),
658  "ERROR: Finite element "
659  << Utility::enum_to_string(base_fe_type.family)
660  << " on geometric element "
661  << Utility::enum_to_string(type)
662  << "\nonly supports FEInterface::max_order = "
663  << FEInterface::max_order(base_fe_type,type)
664  << ", not fe_type.order = "
665  << base_fe_type.order);
666 
667 #ifdef LIBMESH_ENABLE_AMR
668  // Make sure we haven't done more p refinement than we can
669  // handle
670  if (base_fe_type.order + add_p_level*elem->p_level() >
671  FEInterface::max_order(base_fe_type, type))
672  {
673 # ifdef DEBUG
674  libMesh::err << "WARNING: Finite element "
675  << Utility::enum_to_string(base_fe_type.family)
676  << " on geometric element "
677  << Utility::enum_to_string(type) << std::endl
678  << "could not be p refined past FEInterface::max_order = "
679  << FEInterface::max_order(base_fe_type,type)
680  << std::endl;
681 # endif
682  elem->set_p_level(int(FEInterface::max_order(base_fe_type,type))
683  - int(base_fe_type.order));
684  }
685 #endif
686 
687  // Allocate the vertex DOFs
688  for (auto n : elem->node_index_range())
689  {
690  Node & node = elem->node_ref(n);
691 
692  // If we're active on the element then we're active on
693  // its nodes. If we're not then we might *still* be
694  // active on particular constraining nodes.
695  bool active_on_node = active_on_elem;
696  if (!active_on_node)
697  if (auto it = constraining_subdomains.find(&node);
698  it != constraining_subdomains.end())
699  for (auto s : it->second)
700  if (vg_description.active_on_subdomain(s))
701  {
702  active_on_node = true;
703  break;
704  }
705 
706  if (!active_on_node)
707  continue;
708 
709  if (elem->is_vertex(n))
710  {
711  const unsigned int old_node_dofs =
712  node.n_comp_group(sys_num, vg);
713 
714  const unsigned int vertex_dofs =
715  std::max(FEInterface::n_dofs_at_node(fe_type, elem, n, add_p_level),
716  old_node_dofs);
717 
718  // Some discontinuous FEs have no vertex dofs
719  if (vertex_dofs > old_node_dofs)
720  {
721  node.set_n_comp_group(sys_num, vg,
722  vertex_dofs);
723 
724  // Abusing dof_number to set a "this is a
725  // vertex" flag
726  node.set_vg_dof_base(sys_num, vg,
727  vertex_dofs);
728 
729  // libMesh::out << "sys_num,vg,old_node_dofs,vertex_dofs="
730  // << sys_num << ","
731  // << vg << ","
732  // << old_node_dofs << ","
733  // << vertex_dofs << '\n',
734  // node.debug_buffer();
735 
736  // libmesh_assert_equal_to (vertex_dofs, node.n_comp(sys_num, vg));
737  // libmesh_assert_equal_to (vertex_dofs, node.vg_dof_base(sys_num, vg));
738  }
739  }
740  }
741  } // done counting vertex dofs
742 
743  // count edge & face dofs next
744  for (auto & elem : mesh.active_element_ptr_range())
745  {
746  libmesh_assert(elem);
747 
748  // Only number dofs connected to active elements on this
749  // processor and only for variables which are active on on
750  // this element's subdomain or which are active on the
751  // subdomain of a node constrained by this node.
752  const bool active_on_elem =
753  vg_description.active_on_subdomain(elem->subdomain_id());
754 
755  // If there's no way we're active on this element then we're
756  // done
757  if (!active_on_elem && constraining_subdomains_empty)
758  continue;
759 
760  // Allocate the edge and face DOFs
761  for (auto n : elem->node_index_range())
762  {
763  Node & node = elem->node_ref(n);
764 
765  // If we're active on the element then we're active on
766  // its nodes. If we're not then we might *still* be
767  // active on particular constraining nodes.
768  bool active_on_node = active_on_elem;
769  if (!active_on_node)
770  if (auto it = constraining_subdomains.find(&node);
771  it != constraining_subdomains.end())
772  for (auto s : it->second)
773  if (vg_description.active_on_subdomain(s))
774  {
775  active_on_node = true;
776  break;
777  }
778 
779  if (!active_on_node)
780  continue;
781 
782  const unsigned int old_node_dofs =
783  node.n_comp_group(sys_num, vg);
784 
785  const unsigned int vertex_dofs = old_node_dofs?
786  cast_int<unsigned int>(node.vg_dof_base (sys_num,vg)):0;
787 
788  const unsigned int new_node_dofs =
789  FEInterface::n_dofs_at_node(base_fe_type, elem, n, add_p_level);
790 
791  // We've already allocated vertex DOFs
792  if (elem->is_vertex(n))
793  {
794  libmesh_assert_greater_equal (old_node_dofs, vertex_dofs);
795  // //if (vertex_dofs < new_node_dofs)
796  // libMesh::out << "sys_num,vg,old_node_dofs,vertex_dofs,new_node_dofs="
797  // << sys_num << ","
798  // << vg << ","
799  // << old_node_dofs << ","
800  // << vertex_dofs << ","
801  // << new_node_dofs << '\n',
802  // node.debug_buffer();
803 
804  libmesh_assert_greater_equal (vertex_dofs, new_node_dofs);
805  }
806  // We need to allocate the rest
807  else
808  {
809  // If this has no dofs yet, it needs no vertex
810  // dofs, so we just give it edge or face dofs
811  if (!old_node_dofs)
812  {
813  node.set_n_comp_group(sys_num, vg,
814  new_node_dofs);
815  // Abusing dof_number to set a "this has no
816  // vertex dofs" flag
817  if (new_node_dofs)
818  node.set_vg_dof_base(sys_num, vg, 0);
819  }
820 
821  // If this has dofs, but has no vertex dofs,
822  // it may still need more edge or face dofs if
823  // we're p-refined.
824  else if (vertex_dofs == 0)
825  {
826  if (new_node_dofs > old_node_dofs)
827  {
828  node.set_n_comp_group(sys_num, vg,
829  new_node_dofs);
830 
831  node.set_vg_dof_base(sys_num, vg,
832  vertex_dofs);
833  }
834  }
835  // If this is another element's vertex,
836  // add more (non-overlapping) edge/face dofs if
837  // necessary
838  else if (extra_hanging_dofs)
839  {
840  if (new_node_dofs > old_node_dofs - vertex_dofs)
841  {
842  node.set_n_comp_group(sys_num, vg,
843  vertex_dofs + new_node_dofs);
844 
845  node.set_vg_dof_base(sys_num, vg,
846  vertex_dofs);
847  }
848  }
849  // If this is another element's vertex, add any
850  // (overlapping) edge/face dofs if necessary
851  else
852  {
853  libmesh_assert_greater_equal (old_node_dofs, vertex_dofs);
854  if (new_node_dofs > old_node_dofs)
855  {
856  node.set_n_comp_group(sys_num, vg,
857  new_node_dofs);
858 
859  node.set_vg_dof_base (sys_num, vg,
860  vertex_dofs);
861  }
862  }
863  }
864  }
865  // Allocate the element DOFs
866  const unsigned int dofs_per_elem =
867  FEInterface::n_dofs_per_elem(base_fe_type, elem, add_p_level);
868 
869  elem->set_n_comp_group(sys_num, vg, dofs_per_elem);
870 
871  }
872  } // end loop over variable groups
873 
874  // Calling DofMap::reinit() by itself makes little sense,
875  // so we won't bother with nonlocal DofObjects.
876  // Those will be fixed by distribute_dofs
877 
878  //------------------------------------------------------------
879  // Finally, clear all the current DOF indices
880  // (distribute_dofs expects them cleared!)
881  this->invalidate_dofs(mesh);
882 }
883 
884 
885 
887 {
888  const unsigned int sys_num = this->sys_number();
889 
890  // All the nodes
891  for (auto & node : mesh.node_ptr_range())
892  node->invalidate_dofs(sys_num);
893 
894  // All the active elements.
895  for (auto & elem : mesh.active_element_ptr_range())
896  elem->invalidate_dofs(sys_num);
897 }
898 
899 
900 
902 {
904 
905  // we don't want to clear
906  // the coupling matrix!
907  // It should not change...
908  //_dof_coupling->clear();
909  //
910  // But it would be inconsistent to leave our coupling settings
911  // through a clear()...
912  _dof_coupling = nullptr;
913 
914  // Reset ghosting functor statuses
915  {
916  for (const auto & gf : _coupling_functors)
917  {
918  libmesh_assert(gf);
920  }
921  this->_coupling_functors.clear();
922 
923  // Go back to default coupling
924 
925  _default_coupling->set_dof_coupling(this->_dof_coupling);
926  _default_coupling->set_n_levels(this->use_coupled_neighbor_dofs(this->_mesh));
927 
929  }
930 
931 
932  {
933  for (const auto & gf : _algebraic_ghosting_functors)
934  {
935  libmesh_assert(gf);
937  }
938  this->_algebraic_ghosting_functors.clear();
939 
940  // Go back to default send_list generation
941 
942  // _default_evaluating->set_dof_coupling(this->_dof_coupling);
943  _default_evaluating->set_n_levels(1);
945  }
946 
947  this->_shared_functors.clear();
948 
949  _variables.clear();
950  _variable_groups.clear();
951  _var_to_vg.clear();
952  _variable_group_numbers.clear();
953  _first_scalar_df.clear();
954  this->clear_send_list();
955  this->clear_sparsity();
957 
958 #ifdef LIBMESH_ENABLE_AMR
959 
960  _dof_constraints.clear();
961  _stashed_dof_constraints.clear();
964  _n_old_dfs = 0;
965  _first_old_df.clear();
966  _end_old_df.clear();
967  _first_old_scalar_df.clear();
968 
969 #endif
970 
971  _matrices.clear();
972  if (_sc)
973  _sc->clear();
974 }
975 
976 
977 
979 {
980  // This function must be run on all processors at once
981  parallel_object_only();
982 
983  // Log how long it takes to distribute the degrees of freedom
984  LOG_SCOPE("distribute_dofs()", "DofMap");
985 
986  libmesh_assert (mesh.is_prepared());
987 
988  const processor_id_type proc_id = this->processor_id();
989 #ifndef NDEBUG
990  const processor_id_type n_proc = this->n_processors();
991 #endif
992 
993  // libmesh_assert_greater (this->n_variables(), 0);
994  libmesh_assert_less (proc_id, n_proc);
995 
996  // Data structure to ensure we can correctly combine
997  // subdomain-restricted variables with constraining nodes from
998  // different subdomains
999  const std::map<const Node *, std::set<subdomain_id_type>>
1000  constraining_subdomains =
1002 
1003  // re-init in case the mesh has changed
1004  this->reinit(mesh,
1005  constraining_subdomains);
1006 
1007  // By default distribute variables in a
1008  // var-major fashion, but allow run-time
1009  // specification
1010  bool node_major_dofs = libMesh::on_command_line ("--node-major-dofs");
1011 
1012  // The DOF counter, will be incremented as we encounter
1013  // new degrees of freedom
1014  dof_id_type next_free_dof = 0;
1015 
1016  // Clear the send list before we rebuild it
1017  this->clear_send_list();
1018 
1019  // Set temporary DOF indices on this processor
1020  if (node_major_dofs)
1022  (next_free_dof, mesh, constraining_subdomains);
1023  else
1025  (next_free_dof, mesh, constraining_subdomains);
1026 
1027  // Get DOF counts on all processors
1028  const auto n_dofs = this->compute_dof_info(next_free_dof);
1029 
1030  // Clear all the current DOF indices
1031  // (distribute_dofs expects them cleared!)
1032  this->invalidate_dofs(mesh);
1033 
1034  next_free_dof = _first_df[proc_id];
1035 
1036  // Set permanent DOF indices on this processor
1037  if (node_major_dofs)
1039  (next_free_dof, mesh, constraining_subdomains);
1040  else
1042  (next_free_dof, mesh, constraining_subdomains);
1043 
1044  libmesh_assert_equal_to (next_free_dof, _end_df[proc_id]);
1045 
1046  //------------------------------------------------------------
1047  // At this point, all n_comp and dof_number values on local
1048  // DofObjects should be correct, but a DistributedMesh might have
1049  // incorrect values on non-local DofObjects. Let's request the
1050  // correct values from each other processor.
1051 
1052  if (this->n_processors() > 1)
1053  {
1054  this->set_nonlocal_dof_objects(mesh.nodes_begin(),
1055  mesh.nodes_end(),
1057 
1058  this->set_nonlocal_dof_objects(mesh.elements_begin(),
1059  mesh.elements_end(),
1061  }
1062 
1063 #ifdef DEBUG
1064  {
1065  const unsigned int
1066  sys_num = this->sys_number();
1067 
1068  // Processors should all agree on DoF ids for the newly numbered
1069  // system.
1071 
1072  // DoF processor ids should match DofObject processor ids
1073  for (auto & node : mesh.node_ptr_range())
1074  {
1075  DofObject const * const dofobj = node;
1076  const processor_id_type obj_proc_id = dofobj->processor_id();
1077 
1078  for (auto v : make_range(dofobj->n_vars(sys_num)))
1079  for (auto c : make_range(dofobj->n_comp(sys_num,v)))
1080  {
1081  const dof_id_type dofid = dofobj->dof_number(sys_num,v,c);
1082  libmesh_assert_greater_equal (dofid, this->first_dof(obj_proc_id));
1083  libmesh_assert_less (dofid, this->end_dof(obj_proc_id));
1084  }
1085  }
1086 
1087  for (auto & elem : mesh.element_ptr_range())
1088  {
1089  DofObject const * const dofobj = elem;
1090  const processor_id_type obj_proc_id = dofobj->processor_id();
1091 
1092  for (auto v : make_range(dofobj->n_vars(sys_num)))
1093  for (auto c : make_range(dofobj->n_comp(sys_num,v)))
1094  {
1095  const dof_id_type dofid = dofobj->dof_number(sys_num,v,c);
1096  libmesh_assert_greater_equal (dofid, this->first_dof(obj_proc_id));
1097  libmesh_assert_less (dofid, this->end_dof(obj_proc_id));
1098  }
1099  }
1100  }
1101 #endif
1102 
1103  // start finding SCALAR degrees of freedom
1104 #ifdef LIBMESH_ENABLE_AMR
1106 #endif
1107  _first_scalar_df.clear();
1109  dof_id_type current_SCALAR_dof_index = n_dofs - n_SCALAR_dofs();
1110 
1111  // Calculate and cache the initial DoF indices for SCALAR variables.
1112  // This is an O(N_vars) calculation so we want to do it once per
1113  // renumbering rather than once per SCALAR_dof_indices() call
1114 
1115  for (auto v : make_range(this->n_variables()))
1116  if (this->variable(v).type().family == SCALAR)
1117  {
1118  _first_scalar_df[v] = current_SCALAR_dof_index;
1119  current_SCALAR_dof_index += this->variable(v).type().order.get_order();
1120  }
1121 
1122  // Allow our GhostingFunctor objects to reinit if necessary
1123  for (const auto & gf : _algebraic_ghosting_functors)
1124  {
1125  libmesh_assert(gf);
1126  gf->dofmap_reinit();
1127  }
1128 
1129  for (const auto & gf : _coupling_functors)
1130  {
1131  libmesh_assert(gf);
1132  gf->dofmap_reinit();
1133  }
1134 
1135  // Note that in the add_neighbors_to_send_list nodes on processor
1136  // boundaries that are shared by multiple elements are added for
1137  // each element.
1138  this->add_neighbors_to_send_list(mesh);
1139 
1140  // Here we used to clean up that data structure; now System and
1141  // EquationSystems call that for us, after we've added constraint
1142  // dependencies to the send_list too.
1143  // this->sort_send_list ();
1144 
1145  return n_dofs;
1146 }
1147 
1148 
1149 template <typename T, std::enable_if_t<std::is_same_v<T, dof_id_type> ||
1150  std::is_same_v<T, std::vector<dof_id_type>>, int>>
1152  const MeshBase & mesh,
1153  unsigned int var_num) const
1154 {
1155  // Only used if T == dof_id_type to keep track of the greatest dof we've seen
1156  dof_id_type greatest = 0;
1157 
1158  if constexpr (std::is_same_v<T, dof_id_type>)
1159  idx = 0;
1160  else if constexpr (std::is_same_v<T, std::vector<dof_id_type>>)
1161  idx.clear();
1162 
1163  // Count dofs in the *exact* order that distribute_dofs numbered
1164  // them, so that we can assume ascending indices and use push_back
1165  // instead of find+insert.
1166 
1167  const unsigned int sys_num = this->sys_number();
1168 
1169  // If this isn't a SCALAR variable, we need to find all its field
1170  // dofs on the mesh
1171  if (this->variable_type(var_num).family != SCALAR)
1172  {
1173  const Variable & var(this->variable(var_num));
1174 
1175  for (auto & elem : mesh.active_local_element_ptr_range())
1176  {
1177  if (!var.active_on_subdomain(elem->subdomain_id()))
1178  continue;
1179 
1180  // Only count dofs connected to active
1181  // elements on this processor.
1182  const unsigned int n_nodes = elem->n_nodes();
1183 
1184  // First get any new nodal DOFS
1185  for (unsigned int n=0; n<n_nodes; n++)
1186  {
1187  const Node & node = elem->node_ref(n);
1188 
1189  if (node.processor_id() != this->processor_id())
1190  continue;
1191 
1192  const unsigned int n_comp = node.n_comp(sys_num, var_num);
1193  for(unsigned int i=0; i<n_comp; i++)
1194  {
1195  const dof_id_type index = node.dof_number(sys_num,var_num,i);
1196  libmesh_assert (this->local_index(index));
1197 
1198  if constexpr (std::is_same_v<T, dof_id_type>)
1199  {
1200  if (idx == 0 || index > greatest)
1201  { idx++; greatest = index; }
1202  }
1203  else if constexpr (std::is_same_v<T, std::vector<dof_id_type>>)
1204  {
1205  if (idx.empty() || index > idx.back())
1206  idx.push_back(index);
1207  }
1208  }
1209  }
1210 
1211  // Next get any new element DOFS
1212  const unsigned int n_comp = elem->n_comp(sys_num, var_num);
1213  for (unsigned int i=0; i<n_comp; i++)
1214  {
1215  const dof_id_type index = elem->dof_number(sys_num,var_num,i);
1216 
1217  if constexpr (std::is_same_v<T, dof_id_type>)
1218  {
1219  if (idx == 0 || index > greatest)
1220  { idx++; greatest = index; }
1221  }
1222  else if constexpr (std::is_same_v<T, std::vector<dof_id_type>>)
1223  {
1224  if (idx.empty() || index > idx.back())
1225  idx.push_back(index);
1226  }
1227  }
1228  } // done looping over elements
1229 
1230 
1231  // we may have missed assigning DOFs to nodes that we own
1232  // but to which we have no connected elements matching our
1233  // variable restriction criterion. this will happen, for example,
1234  // if variable V is restricted to subdomain S. We may not own
1235  // any elements which live in S, but we may own nodes which are
1236  // *connected* to elements which do. in this scenario these nodes
1237  // will presently have unnumbered DOFs. we need to take care of
1238  // them here since we own them and no other processor will touch them.
1239  for (const auto & node : mesh.local_node_ptr_range())
1240  {
1241  libmesh_assert(node);
1242 
1243  const unsigned int n_comp = node->n_comp(sys_num, var_num);
1244  for (unsigned int i=0; i<n_comp; i++)
1245  {
1246  const dof_id_type index = node->dof_number(sys_num,var_num,i);
1247 
1248  if constexpr (std::is_same_v<T, dof_id_type>)
1249  {
1250  if (idx == 0 || index > greatest)
1251  { idx++; greatest = index; }
1252  }
1253  else if constexpr (std::is_same_v<T, std::vector<dof_id_type>>)
1254  {
1255  if (idx.empty() || index > idx.back())
1256  idx.push_back(index);
1257  }
1258  }
1259  }
1260  }
1261  // Otherwise, count up the SCALAR dofs, if we're on the processor
1262  // that holds this SCALAR variable
1263  else if (this->processor_id() == (this->n_processors()-1))
1264  {
1265  std::vector<dof_id_type> di_scalar;
1266  this->SCALAR_dof_indices(di_scalar,var_num);
1267 
1268  if constexpr (std::is_same_v<T, dof_id_type>)
1269  idx += std::distance(di_scalar.begin(), di_scalar.end());
1270  else if constexpr (std::is_same_v<T, std::vector<dof_id_type>>)
1271  idx.insert(idx.end(), di_scalar.begin(), di_scalar.end());
1272  }
1273 }
1274 
1276  const MeshBase &,
1277  unsigned int) const;
1278 
1279 template void DofMap::local_variable_indices(std::vector<dof_id_type> &,
1280  const MeshBase &,
1281  unsigned int) const;
1282 
1283 
1284 std::map<const Node *, std::set<subdomain_id_type>>
1286 {
1287  std::map<const Node *, std::set<subdomain_id_type>> constraining_subdomains;
1288  const auto & constraint_rows = _mesh.get_constraint_rows();
1289 
1290  // We can't just loop over constraint rows here because we need
1291  // element subdomain ids for the constrained nodes, but we don't
1292  // want an extra loop if there are no constraint rows.
1293  if (!constraint_rows.empty())
1294  for (auto & elem : _mesh.active_element_ptr_range())
1295  {
1296  const subdomain_id_type sbdid = elem->subdomain_id();
1297 
1298  for (const Node & node : elem->node_ref_range())
1299  {
1300  if (auto it = constraint_rows.find(&node);
1301  it != constraint_rows.end())
1302  {
1303  for (const auto & [pr, val] : it->second)
1304  {
1305  const Node * spline_node =
1306  pr.first->node_ptr(pr.second);
1307 
1308  constraining_subdomains[spline_node].insert(sbdid);
1309  }
1310  }
1311  }
1312  }
1313 
1314  return constraining_subdomains;
1315 }
1316 
1317 
1319  (dof_id_type & next_free_dof,
1320  MeshBase & mesh,
1321  const std::map<const Node *, std::set<subdomain_id_type>> &
1322  constraining_subdomains)
1323 {
1324  const unsigned int sys_num = this->sys_number();
1325  const unsigned int n_var_groups = this->n_variable_groups();
1326 
1327  // This is the common case and we want to optimize for it
1328  const bool constraining_subdomains_empty =
1329  constraining_subdomains.empty();
1330 
1331  // Our numbering here must be kept consistent with the numbering
1332  // scheme assumed by DofMap::local_variable_indices!
1333 
1334  //-------------------------------------------------------------------------
1335  // First count and assign temporary numbers to local dofs
1336  for (auto & elem : mesh.active_local_element_ptr_range())
1337  {
1338  // Only number dofs connected to active
1339  // elements on this processor.
1340  const unsigned int n_nodes = elem->n_nodes();
1341 
1342  const subdomain_id_type sbdid = elem->subdomain_id();
1343 
1344  // First number the nodal DOFS
1345  for (unsigned int n=0; n<n_nodes; n++)
1346  {
1347  Node & node = elem->node_ref(n);
1348 
1349  for (unsigned vg=0; vg<n_var_groups; vg++)
1350  {
1351  const VariableGroup & vg_description(this->variable_group(vg));
1352 
1353  if (vg_description.type().family == SCALAR)
1354  continue;
1355 
1356  bool active_on_node =
1357  vg_description.active_on_subdomain(sbdid);
1358 
1359  // Are we at least active indirectly here?
1360  if (!active_on_node && !constraining_subdomains_empty)
1361  if (auto it = constraining_subdomains.find(&node);
1362  it != constraining_subdomains.end())
1363  for (auto s : it->second)
1364  if (vg_description.active_on_subdomain(s))
1365  {
1366  active_on_node = true;
1367  break;
1368  }
1369 
1370  if (active_on_node)
1371  {
1372  // assign dof numbers (all at once) if this is
1373  // our node and if they aren't already there
1374  if ((node.n_comp_group(sys_num,vg) > 0) &&
1375  (node.processor_id() == this->processor_id()) &&
1376  (node.vg_dof_base(sys_num,vg) ==
1378  {
1379  node.set_vg_dof_base(sys_num, vg,
1380  next_free_dof);
1381  next_free_dof += (vg_description.n_variables()*
1382  node.n_comp_group(sys_num,vg));
1383  //node.debug_buffer();
1384  }
1385  }
1386  }
1387  }
1388 
1389  // Now number the element DOFS
1390  for (unsigned vg=0; vg<n_var_groups; vg++)
1391  {
1392  const VariableGroup & vg_description(this->variable_group(vg));
1393 
1394  if ((vg_description.type().family != SCALAR) &&
1395  (vg_description.active_on_subdomain(elem->subdomain_id())))
1396  if (elem->n_comp_group(sys_num,vg) > 0)
1397  {
1398  libmesh_assert_equal_to (elem->vg_dof_base(sys_num,vg),
1400 
1401  elem->set_vg_dof_base(sys_num,
1402  vg,
1403  next_free_dof);
1404 
1405  next_free_dof += (vg_description.n_variables()*
1406  elem->n_comp_group(sys_num,vg));
1407  }
1408  }
1409  } // done looping over elements
1410 
1411 
1412  // we may have missed assigning DOFs to nodes that we own
1413  // but to which we have no connected elements matching our
1414  // variable restriction criterion. this will happen, for example,
1415  // if variable V is restricted to subdomain S. We may not own
1416  // any elements which live in S, but we may own nodes which are
1417  // *connected* to elements which do. in this scenario these nodes
1418  // will presently have unnumbered DOFs. we need to take care of
1419  // them here since we own them and no other processor will touch them.
1420  for (auto & node : mesh.local_node_ptr_range())
1421  for (unsigned vg=0; vg<n_var_groups; vg++)
1422  {
1423  const VariableGroup & vg_description(this->variable_group(vg));
1424 
1425  if (node->n_comp_group(sys_num,vg))
1426  if (node->vg_dof_base(sys_num,vg) == DofObject::invalid_id)
1427  {
1428  node->set_vg_dof_base (sys_num,
1429  vg,
1430  next_free_dof);
1431 
1432  next_free_dof += (vg_description.n_variables()*
1433  node->n_comp(sys_num,vg));
1434  }
1435  }
1436 
1437  this->distribute_scalar_dofs(next_free_dof);
1438 
1439 #ifdef DEBUG
1440  this->assert_no_nodes_missed(mesh);
1441 #endif // DEBUG
1442 }
1443 
1444 
1445 
1447  (dof_id_type & next_free_dof,
1448  MeshBase & mesh,
1449  const std::map<const Node *, std::set<subdomain_id_type>> &
1450  constraining_subdomains)
1451 {
1452  const unsigned int sys_num = this->sys_number();
1453  const unsigned int n_var_groups = this->n_variable_groups();
1454 
1455  // This is the common case and we want to optimize for it
1456  const bool constraining_subdomains_empty =
1457  constraining_subdomains.empty();
1458 
1459  // Our numbering here must be kept consistent with the numbering
1460  // scheme assumed by DofMap::local_variable_indices!
1461 
1462  //-------------------------------------------------------------------------
1463  // First count and assign temporary numbers to local dofs
1464  for (unsigned vg=0; vg<n_var_groups; vg++)
1465  {
1466  const VariableGroup & vg_description(this->variable_group(vg));
1467 
1468  const unsigned int n_vars_in_group = vg_description.n_variables();
1469 
1470  // Skip the SCALAR dofs
1471  if (vg_description.type().family == SCALAR)
1472  continue;
1473 
1474  for (auto & elem : mesh.active_local_element_ptr_range())
1475  {
1476  // Only number dofs connected to active elements on this
1477  // processor and only for variables which are active on on
1478  // this element's subdomain or which are active on the
1479  // subdomain of a node constrained by this node.
1480  const bool active_on_elem =
1481  vg_description.active_on_subdomain(elem->subdomain_id());
1482 
1483  // If there's no way we're active on this element then we're
1484  // done
1485  if (!active_on_elem && constraining_subdomains_empty)
1486  continue;
1487 
1488  const unsigned int n_nodes = elem->n_nodes();
1489 
1490  // First number the nodal DOFS
1491  for (unsigned int n=0; n<n_nodes; n++)
1492  {
1493  Node & node = elem->node_ref(n);
1494 
1495  bool active_on_node = active_on_elem;
1496  if (!active_on_node)
1497  if (auto it = constraining_subdomains.find(&node);
1498  it != constraining_subdomains.end())
1499  for (auto s : it->second)
1500  if (vg_description.active_on_subdomain(s))
1501  {
1502  active_on_node = true;
1503  break;
1504  }
1505 
1506  if (!active_on_node)
1507  continue;
1508 
1509  // assign dof numbers (all at once) if this is
1510  // our node and if they aren't already there
1511  if ((node.n_comp_group(sys_num,vg) > 0) &&
1512  (node.processor_id() == this->processor_id()) &&
1513  (node.vg_dof_base(sys_num,vg) ==
1515  {
1516  node.set_vg_dof_base(sys_num, vg, next_free_dof);
1517 
1518  next_free_dof += (n_vars_in_group*
1519  node.n_comp_group(sys_num,vg));
1520  }
1521  }
1522 
1523  // Now number the element DOFS
1524  if (elem->n_comp_group(sys_num,vg) > 0)
1525  {
1526  libmesh_assert_equal_to (elem->vg_dof_base(sys_num,vg),
1528 
1529  elem->set_vg_dof_base(sys_num,
1530  vg,
1531  next_free_dof);
1532 
1533  next_free_dof += (n_vars_in_group*
1534  elem->n_comp_group(sys_num,vg));
1535  }
1536  } // end loop on elements
1537 
1538  // we may have missed assigning DOFs to nodes that we own
1539  // but to which we have no connected elements matching our
1540  // variable restriction criterion. this will happen, for example,
1541  // if variable V is restricted to subdomain S. We may not own
1542  // any elements which live in S, but we may own nodes which are
1543  // *connected* to elements which do. in this scenario these nodes
1544  // will presently have unnumbered DOFs. we need to take care of
1545  // them here since we own them and no other processor will touch them.
1546  for (auto & node : mesh.local_node_ptr_range())
1547  if (node->n_comp_group(sys_num,vg))
1548  if (node->vg_dof_base(sys_num,vg) == DofObject::invalid_id)
1549  {
1550  node->set_vg_dof_base (sys_num,
1551  vg,
1552  next_free_dof);
1553 
1554  next_free_dof += (n_vars_in_group*
1555  node->n_comp_group(sys_num,vg));
1556  }
1557  } // end loop on variable groups
1558 
1559  this->distribute_scalar_dofs(next_free_dof);
1560 
1561 #ifdef DEBUG
1562  this->assert_no_nodes_missed(mesh);
1563 #endif
1564 }
1565 
1566 
1567 
1569 {
1570  this->_n_SCALAR_dofs = 0;
1571  for (auto vg : make_range(this->n_variable_groups()))
1572  {
1573  const VariableGroup & vg_description(this->variable_group(vg));
1574 
1575  if (vg_description.type().family == SCALAR)
1576  {
1577  this->_n_SCALAR_dofs += (vg_description.n_variables()*
1578  vg_description.type().order.get_order());
1579  continue;
1580  }
1581  }
1582 
1583  // Only increment next_free_dof if we're on the processor
1584  // that holds this SCALAR variable
1585  if (this->processor_id() == (this->n_processors()-1))
1586  next_free_dof += _n_SCALAR_dofs;
1587 }
1588 
1589 
1590 
1591 #ifdef DEBUG
1593 {
1594  MeshTools::libmesh_assert_valid_procids<Node>(mesh);
1595 
1596  for (auto & node : mesh.local_node_ptr_range())
1597  {
1598  unsigned int n_var_g = node->n_var_groups(this->sys_number());
1599  for (unsigned int vg=0; vg != n_var_g; ++vg)
1600  {
1601  unsigned int n_comp_g =
1602  node->n_comp_group(this->sys_number(), vg);
1603  dof_id_type my_first_dof = n_comp_g ?
1604  node->vg_dof_base(this->sys_number(), vg) : 0;
1605  libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
1606  }
1607  }
1608 }
1609 #endif // DEBUG
1610 
1611 
1612 void
1613 DofMap::
1615  CouplingMatricesSet & temporary_coupling_matrices,
1616  const GhostingFunctorIterator & gf_begin,
1617  const GhostingFunctorIterator & gf_end,
1618  const MeshBase::const_element_iterator & elems_begin,
1619  const MeshBase::const_element_iterator & elems_end,
1621 {
1622  for (const auto & gf : as_range(gf_begin, gf_end))
1623  {
1624  GhostingFunctor::map_type more_elements_to_ghost;
1625 
1626  libmesh_assert(gf);
1627  (*gf)(elems_begin, elems_end, p, more_elements_to_ghost);
1628 
1629  // A GhostingFunctor should only return active elements, but
1630  // I forgot to *document* that, so let's go as easy as we
1631  // can on functors that return inactive elements.
1632 #if defined(LIBMESH_ENABLE_DEPRECATED) && defined(LIBMESH_ENABLE_AMR)
1633  std::vector<std::pair<const Elem*, const CouplingMatrix*>> children_to_couple;
1634  for (auto it = more_elements_to_ghost.begin();
1635  it != more_elements_to_ghost.end();)
1636  {
1637  const Elem * elem = it->first;
1638  if (!elem->active())
1639  {
1640  libmesh_deprecated();
1641  std::vector<const Elem*> children_to_ghost;
1642  elem->active_family_tree(children_to_ghost,
1643  /*reset=*/ false);
1644  for (const Elem * child : children_to_ghost)
1645  if (child->processor_id() != p)
1646  children_to_couple.emplace_back(child, it->second);
1647 
1648  it = more_elements_to_ghost.erase(it);
1649  }
1650  else
1651  ++it;
1652  }
1653  more_elements_to_ghost.insert(children_to_couple.begin(),
1654  children_to_couple.end());
1655 #endif
1656 
1657  for (const auto & [elem, elem_cm] : more_elements_to_ghost)
1658  {
1659  // At this point we should only have active elements, even
1660  // if we had to fix up gf output to get here.
1661  libmesh_assert(elem->active());
1662 
1663  if (const auto existing_it = elements_to_ghost.find(elem);
1664  existing_it == elements_to_ghost.end())
1665  elements_to_ghost.emplace(elem, elem_cm);
1666  else
1667  {
1668  if (existing_it->second)
1669  {
1670  if (elem_cm)
1671  {
1672  // If this isn't already a temporary
1673  // then we need to make one so we'll
1674  // have a non-const matrix to merge
1675  if (temporary_coupling_matrices.empty() ||
1676  !temporary_coupling_matrices.count(existing_it->second))
1677  {
1678  // Make copy. This just calls the
1679  // compiler-generated copy constructor
1680  // because the CouplingMatrix class does not
1681  // define a custom copy constructor.
1682  auto result_pr = temporary_coupling_matrices.insert(std::make_unique<CouplingMatrix>(*existing_it->second));
1683  existing_it->second = result_pr.first->get();
1684  }
1685 
1686  // Merge elem_cm into existing CouplingMatrix
1687  const_cast<CouplingMatrix &>(*existing_it->second) &= *elem_cm;
1688  }
1689  else // elem_cm == nullptr
1690  {
1691  // Any existing_it matrix merged with a full
1692  // matrix (symbolized as nullptr) gives another
1693  // full matrix (symbolizable as nullptr).
1694 
1695  // So if existing_it->second is a temporary then
1696  // we don't need it anymore; we might as well
1697  // remove it to keep the set of temporaries
1698  // small.
1699  if (const auto temp_it = temporary_coupling_matrices.find(existing_it->second);
1700  temp_it != temporary_coupling_matrices.end())
1701  temporary_coupling_matrices.erase(temp_it);
1702 
1703  existing_it->second = nullptr;
1704  }
1705  }
1706  // else we have a nullptr already, then we have a full
1707  // coupling matrix, already, and merging with anything
1708  // else won't change that, so we're done.
1709  }
1710  }
1711  }
1712 }
1713 
1714 
1715 
1717 {
1718  LOG_SCOPE("add_neighbors_to_send_list()", "DofMap");
1719 
1720  // Return immediately if there's no ghost data
1721  if (this->n_processors() == 1)
1722  return;
1723 
1724  const unsigned int n_var = this->n_variables();
1725 
1726  MeshBase::const_element_iterator local_elem_it
1727  = mesh.active_local_elements_begin();
1728  const MeshBase::const_element_iterator local_elem_end
1729  = mesh.active_local_elements_end();
1730 
1731  GhostingFunctor::map_type elements_to_send;
1732  DofMap::CouplingMatricesSet temporary_coupling_matrices;
1733 
1734  // We need to add dofs to the send list if they've been directly
1735  // requested by an algebraic ghosting functor or they've been
1736  // indirectly requested by a coupling functor.
1737  this->merge_ghost_functor_outputs(elements_to_send,
1738  temporary_coupling_matrices,
1741  local_elem_it, local_elem_end, mesh.processor_id());
1742 
1743  this->merge_ghost_functor_outputs(elements_to_send,
1744  temporary_coupling_matrices,
1745  this->coupling_functors_begin(),
1746  this->coupling_functors_end(),
1747  local_elem_it, local_elem_end, mesh.processor_id());
1748 
1749  // Making a list of non-zero coupling matrix columns is an
1750  // O(N_var^2) operation. We cache it so we only have to do it once
1751  // per CouplingMatrix and not once per element.
1752  std::map<const CouplingMatrix *, std::vector<unsigned int>>
1753  column_variable_lists;
1754 
1755  for (const auto & [partner, ghost_coupling] : elements_to_send)
1756  {
1757  // We asked ghosting functors not to give us local elements
1758  libmesh_assert_not_equal_to
1759  (partner->processor_id(), this->processor_id());
1760 
1761  // Loop over any present coupling matrix column variables if we
1762  // have a coupling matrix, or just add all variables to
1763  // send_list if not.
1764  if (ghost_coupling)
1765  {
1766  libmesh_assert_equal_to (ghost_coupling->size(), n_var);
1767 
1768  // Try to find a cached list of column variables.
1769  std::map<const CouplingMatrix *, std::vector<unsigned int>>::const_iterator
1770  column_variable_list = column_variable_lists.find(ghost_coupling);
1771 
1772  // If we didn't find it, then we need to create it.
1773  if (column_variable_list == column_variable_lists.end())
1774  {
1775  auto inserted_variable_list_pair =
1776  column_variable_lists.emplace(ghost_coupling, std::vector<unsigned int>());
1777  column_variable_list = inserted_variable_list_pair.first;
1778 
1779  std::vector<unsigned int> & new_variable_list =
1780  inserted_variable_list_pair.first->second;
1781 
1782  std::vector<unsigned char> has_variable(n_var, false);
1783 
1784  for (unsigned int vi = 0; vi != n_var; ++vi)
1785  {
1786  ConstCouplingRow ccr(vi, *ghost_coupling);
1787 
1788  for (const auto & vj : ccr)
1789  has_variable[vj] = true;
1790  }
1791  for (unsigned int vj = 0; vj != n_var; ++vj)
1792  {
1793  if (has_variable[vj])
1794  new_variable_list.push_back(vj);
1795  }
1796  }
1797 
1798  const std::vector<unsigned int> & variable_list =
1799  column_variable_list->second;
1800 
1801  for (const auto & vj : variable_list)
1802  {
1803  std::vector<dof_id_type> di;
1804  this->dof_indices (partner, di, vj);
1805 
1806  // Insert the remote DOF indices into the send list
1807  for (auto d : di)
1808  if (d != DofObject::invalid_id &&
1809  !this->local_index(d))
1810  {
1811  libmesh_assert_less(d, this->n_dofs());
1812  _send_list.push_back(d);
1813  }
1814  }
1815  }
1816  else
1817  {
1818  std::vector<dof_id_type> di;
1819  this->dof_indices (partner, di);
1820 
1821  // Insert the remote DOF indices into the send list
1822  for (const auto & dof : di)
1823  if (dof != DofObject::invalid_id &&
1824  !this->local_index(dof))
1825  {
1826  libmesh_assert_less(dof, this->n_dofs());
1827  _send_list.push_back(dof);
1828  }
1829  }
1830 
1831  }
1832 
1833  // We're now done with any merged coupling matrices we had to create.
1834  temporary_coupling_matrices.clear();
1835 
1836  //-------------------------------------------------------------------------
1837  // Our coupling functors added dofs from neighboring elements to the
1838  // send list, but we may still need to add non-local dofs from local
1839  // elements.
1840  //-------------------------------------------------------------------------
1841 
1842  // Loop over the active local elements, adding all active elements
1843  // that neighbor an active local element to the send list.
1844  for ( ; local_elem_it != local_elem_end; ++local_elem_it)
1845  {
1846  const Elem * elem = *local_elem_it;
1847 
1848  std::vector<dof_id_type> di;
1849  this->dof_indices (elem, di);
1850 
1851  // Insert the remote DOF indices into the send list
1852  for (const auto & dof : di)
1853  if (dof != DofObject::invalid_id &&
1854  !this->local_index(dof))
1855  {
1856  libmesh_assert_less(dof, this->n_dofs());
1857  _send_list.push_back(dof);
1858  }
1859  }
1860 }
1861 
1862 
1863 
1865 {
1866  LOG_SCOPE("prepare_send_list()", "DofMap");
1867 
1868  // Return immediately if there's no ghost data
1869  if (this->n_processors() == 1)
1870  return;
1871 
1872  // Check to see if we have any extra stuff to add to the send_list
1874  {
1875  if (_augment_send_list)
1876  {
1877  libmesh_here();
1878  libMesh::out << "WARNING: You have specified both an extra send list function and object.\n"
1879  << " Are you sure this is what you meant to do??"
1880  << std::endl;
1881  }
1882 
1884  }
1885 
1886  if (_augment_send_list)
1888 
1889  // First sort the send list. After this
1890  // duplicated elements will be adjacent in the
1891  // vector
1892  std::sort(_send_list.begin(), _send_list.end());
1893 
1894  // Now use std::unique to remove duplicate entries
1895  std::vector<dof_id_type>::iterator new_end =
1896  std::unique (_send_list.begin(), _send_list.end());
1897 
1898  // Remove the end of the send_list. Use the "swap trick"
1899  // from Effective STL
1900  std::vector<dof_id_type> (_send_list.begin(), new_end).swap (_send_list);
1901 
1902  // Make sure the send list has nothing invalid in it.
1903  libmesh_assert(_send_list.empty() || _send_list.back() < this->n_dofs());
1904 }
1905 
1907 {
1908  this->clear_send_list();
1909  this->add_neighbors_to_send_list(mesh);
1910 
1911 #ifdef LIBMESH_ENABLE_CONSTRAINTS
1912  // This is assuming that we only need to recommunicate
1913  // the constraints and no new ones have been added since
1914  // a previous call to reinit_constraints.
1915  this->process_constraints(mesh);
1916 #endif
1917  this->prepare_send_list();
1918 }
1919 
1920 void DofMap::set_implicit_neighbor_dofs(bool implicit_neighbor_dofs)
1921 {
1923  _implicit_neighbor_dofs = implicit_neighbor_dofs;
1924 }
1925 
1927 {
1929 }
1930 
1931 
1932 bool DofMap::use_coupled_neighbor_dofs(const MeshBase & /*mesh*/) const
1933 {
1934  // If we were asked on the command line, then we need to
1935  // include sensitivities between neighbor degrees of freedom
1936  bool implicit_neighbor_dofs =
1937  libMesh::on_command_line ("--implicit-neighbor-dofs");
1938 
1939  // If the user specifies --implicit-neighbor-dofs 0, then
1940  // presumably he knows what he is doing and we won't try to
1941  // automatically turn it on even when all the variables are
1942  // discontinuous.
1943  if (implicit_neighbor_dofs)
1944  {
1945  // No flag provided defaults to 'true'
1946  int flag = 1;
1947  flag = libMesh::command_line_next ("--implicit-neighbor-dofs", flag);
1948 
1949  if (!flag)
1950  {
1951  // The user said --implicit-neighbor-dofs 0, so he knows
1952  // what he is doing and really doesn't want it.
1953  return false;
1954  }
1955  }
1956 
1957  // Possibly override the commandline option, if set_implicit_neighbor_dofs
1958  // has been called.
1960  {
1961  implicit_neighbor_dofs = _implicit_neighbor_dofs;
1962 
1963  // Again, if the user explicitly says implicit_neighbor_dofs = false,
1964  // then we return here.
1965  if (!implicit_neighbor_dofs)
1966  return false;
1967  }
1968 
1969  // Look at all the variables in this system. If every one is
1970  // discontinuous then the user must be doing DG/FVM, so be nice
1971  // and force implicit_neighbor_dofs=true.
1972  {
1973  bool all_discontinuous_dofs = true;
1974 
1975  for (auto var : make_range(this->n_variables()))
1977  all_discontinuous_dofs = false;
1978 
1979  if (all_discontinuous_dofs)
1980  implicit_neighbor_dofs = true;
1981  }
1982 
1983  return implicit_neighbor_dofs;
1984 }
1985 
1986 
1987 
1989 {
1991 
1992  // It is possible that some \p SparseMatrix implementations want to
1993  // see the sparsity pattern before we throw it away. If so, we
1994  // share a view of its arrays, and we pass it in to the matrices.
1995  for (const auto & mat : _matrices)
1996  {
1997  mat->attach_sparsity_pattern (*_sp);
1999  mat->update_sparsity_pattern (_sp->get_sparsity_pattern());
2000  }
2001  // If we don't need the full sparsity pattern anymore, free the
2002  // parts of it we don't need.
2004  _sp->clear_full_sparsity();
2005 }
2006 
2007 
2008 
2010 {
2011  _sp.reset();
2012 }
2013 
2014 
2015 
2017 {
2020 }
2021 
2022 
2023 
2025 {
2026  this->add_coupling_functor(this->default_coupling());
2028 }
2029 
2030 
2031 
2032 void
2034  bool to_mesh)
2035 {
2036  // We used to implicitly support duplicate inserts to std::set
2037 #ifdef LIBMESH_ENABLE_DEPRECATED
2038  _coupling_functors.erase
2039  (std::remove(_coupling_functors.begin(),
2040  _coupling_functors.end(),
2041  &coupling_functor),
2042  _coupling_functors.end());
2043 #endif
2044 
2045  // We shouldn't have two copies of the same functor
2046  libmesh_assert(std::find(_coupling_functors.begin(),
2047  _coupling_functors.end(),
2048  &coupling_functor) ==
2049  _coupling_functors.end());
2050 
2051  _coupling_functors.push_back(&coupling_functor);
2052  coupling_functor.set_mesh(&_mesh);
2053  if (to_mesh)
2054  _mesh.add_ghosting_functor(coupling_functor);
2055 }
2056 
2057 
2058 
2059 void
2061 {
2062  auto raw_it = std::find(_coupling_functors.begin(),
2063  _coupling_functors.end(), &coupling_functor);
2064 
2065 #ifndef LIBMESH_ENABLE_DEPRECATED
2066  // We shouldn't be trying to remove a functor that isn't there
2067  libmesh_assert(raw_it != _coupling_functors.end());
2068 #else
2069  // Our old API supported trying to remove a functor that isn't there
2070  if (raw_it != _coupling_functors.end())
2071 #endif
2072  _coupling_functors.erase(raw_it);
2073 
2074  // We shouldn't have had two copies of the same functor
2075  libmesh_assert(std::find(_coupling_functors.begin(),
2076  _coupling_functors.end(),
2077  &coupling_functor) ==
2078  _coupling_functors.end());
2079 
2080  _mesh.remove_ghosting_functor(coupling_functor);
2081 
2082  if (const auto it = _shared_functors.find(&coupling_functor);
2083  it != _shared_functors.end())
2084  _shared_functors.erase(it);
2085 }
2086 
2087 
2088 
2089 void
2091  bool to_mesh)
2092 {
2093  // We used to implicitly support duplicate inserts to std::set
2094 #ifdef LIBMESH_ENABLE_DEPRECATED
2096  (std::remove(_algebraic_ghosting_functors.begin(),
2098  &evaluable_functor),
2100 #endif
2101 
2102  // We shouldn't have two copies of the same functor
2103  libmesh_assert(std::find(_algebraic_ghosting_functors.begin(),
2105  &evaluable_functor) ==
2107 
2108  _algebraic_ghosting_functors.push_back(&evaluable_functor);
2109  evaluable_functor.set_mesh(&_mesh);
2110  if (to_mesh)
2111  _mesh.add_ghosting_functor(evaluable_functor);
2112 }
2113 
2114 
2115 
2116 void
2118 {
2119  auto raw_it = std::find(_algebraic_ghosting_functors.begin(),
2121  &evaluable_functor);
2122 
2123 #ifndef LIBMESH_ENABLE_DEPRECATED
2124  // We shouldn't be trying to remove a functor that isn't there
2126 #else
2127  // Our old API supported trying to remove a functor that isn't there
2128  if (raw_it != _algebraic_ghosting_functors.end())
2129 #endif
2130  _algebraic_ghosting_functors.erase(raw_it);
2131 
2132  // We shouldn't have had two copies of the same functor
2133  libmesh_assert(std::find(_algebraic_ghosting_functors.begin(),
2135  &evaluable_functor) ==
2137 
2138  _mesh.remove_ghosting_functor(evaluable_functor);
2139 
2140  if (const auto it = _shared_functors.find(&evaluable_functor);
2141  it != _shared_functors.end())
2142  _shared_functors.erase(it);
2143 }
2144 
2145 
2146 
2148  const std::vector<dof_id_type> & dof_indices_in,
2149  DenseVectorBase<Number> & Ue) const
2150 {
2151  const unsigned int n_original_dofs = dof_indices_in.size();
2152 
2153 #ifdef LIBMESH_ENABLE_AMR
2154 
2155  // Trivial mapping
2156  libmesh_assert_equal_to (dof_indices_in.size(), Ue.size());
2157  bool has_constrained_dofs = false;
2158 
2159  for (unsigned int il=0; il != n_original_dofs; ++il)
2160  {
2161  const dof_id_type ig = dof_indices_in[il];
2162 
2163  if (this->is_constrained_dof (ig)) has_constrained_dofs = true;
2164 
2165  libmesh_assert_less (ig, Ug.size());
2166 
2167  Ue.el(il) = Ug(ig);
2168  }
2169 
2170  // If the element has any constrained DOFs then we need
2171  // to account for them in the mapping. This will handle
2172  // the case that the input vector is not constrained.
2173  if (has_constrained_dofs)
2174  {
2175  // Copy the input DOF indices.
2176  std::vector<dof_id_type> constrained_dof_indices(dof_indices_in);
2177 
2180 
2181  this->build_constraint_matrix_and_vector (C, H, constrained_dof_indices);
2182 
2183  libmesh_assert_equal_to (dof_indices_in.size(), C.m());
2184  libmesh_assert_equal_to (constrained_dof_indices.size(), C.n());
2185 
2186  // zero-out Ue
2187  Ue.zero();
2188 
2189  // compute Ue = C Ug, with proper mapping.
2190  for (unsigned int i=0; i != n_original_dofs; i++)
2191  {
2192  Ue.el(i) = H(i);
2193 
2194  const unsigned int n_constrained =
2195  cast_int<unsigned int>(constrained_dof_indices.size());
2196  for (unsigned int j=0; j<n_constrained; j++)
2197  {
2198  const dof_id_type jg = constrained_dof_indices[j];
2199 
2200  // If Ug is a serial or ghosted vector, then this assert is
2201  // overzealous. If Ug is a parallel vector, then this assert
2202  // is redundant.
2203  // libmesh_assert ((jg >= Ug.first_local_index()) &&
2204  // (jg < Ug.last_local_index()));
2205 
2206  Ue.el(i) += C(i,j)*Ug(jg);
2207  }
2208  }
2209  }
2210 
2211 #else
2212 
2213  // Trivial mapping
2214 
2215  libmesh_assert_equal_to (n_original_dofs, Ue.size());
2216 
2217  for (unsigned int il=0; il<n_original_dofs; il++)
2218  {
2219  const dof_id_type ig = dof_indices_in[il];
2220 
2221  libmesh_assert ((ig >= Ug.first_local_index()) && (ig < Ug.last_local_index()));
2222 
2223  Ue.el(il) = Ug(ig);
2224  }
2225 
2226 #endif
2227 }
2228 
2229 void DofMap::dof_indices (const Elem * const elem,
2230  std::vector<dof_id_type> & di) const
2231 {
2232  // We now allow elem==nullptr to request just SCALAR dofs
2233  // libmesh_assert(elem);
2234 
2235  // If we are asking for current indices on an element, it ought to
2236  // be an active element (or a temporary side, which also thinks it's
2237  // active)
2238  libmesh_assert(!elem || elem->active());
2239 
2240  // dof_indices() is a relatively light-weight function that is
2241  // called millions of times in normal codes. Therefore, it is not a
2242  // good candidate for logging, since the cost of the logging code
2243  // itself is roughly on par with the time required to call
2244  // dof_indices().
2245  // LOG_SCOPE("dof_indices()", "DofMap");
2246 
2247  // Clear the DOF indices vector
2248  di.clear();
2249 
2250  const unsigned int n_var_groups = this->n_variable_groups();
2251 
2252 #ifdef DEBUG
2253  // Check that sizes match in DEBUG mode
2254  std::size_t tot_size = 0;
2255 #endif
2256 
2257  if (elem && elem->type() == TRI3SUBDIVISION)
2258  {
2259  // Subdivision surface FE require the 1-ring around elem
2260  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
2261 
2262  // Ghost subdivision elements have no real dofs
2263  if (!sd_elem->is_ghost())
2264  {
2265  // Determine the nodes contributing to element elem
2266  std::vector<const Node *> elem_nodes;
2267  MeshTools::Subdivision::find_one_ring(sd_elem, elem_nodes);
2268 
2269  // Get the dof numbers
2270  for (unsigned int vg=0; vg<n_var_groups; vg++)
2271  {
2272  const VariableGroup & var = this->variable_group(vg);
2273  const unsigned int vars_in_group = var.n_variables();
2274 
2275  if (var.type().family == SCALAR &&
2276  var.active_on_subdomain(elem->subdomain_id()))
2277  {
2278  for (unsigned int vig=0; vig != vars_in_group; ++vig)
2279  {
2280 #ifdef DEBUG
2281  tot_size += var.type().order;
2282 #endif
2283  std::vector<dof_id_type> di_new;
2284  this->SCALAR_dof_indices(di_new,var.number(vig));
2285  di.insert( di.end(), di_new.begin(), di_new.end());
2286  }
2287  }
2288  else
2289  for (unsigned int vig=0; vig != vars_in_group; ++vig)
2290  {
2291  _dof_indices(*elem, elem->p_level(), di, vg, vig,
2292  elem_nodes.data(),
2293  cast_int<unsigned int>(elem_nodes.size()),
2294  var.number(vig)
2295 #ifdef DEBUG
2296  , tot_size
2297 #endif
2298  );
2299  }
2300  }
2301  }
2302 
2303  return;
2304  }
2305 
2306  // Get the dof numbers for each variable
2307  const unsigned int n_nodes = elem ? elem->n_nodes() : 0;
2308  for (unsigned int vg=0; vg<n_var_groups; vg++)
2309  {
2310  const VariableGroup & var = this->variable_group(vg);
2311  const unsigned int vars_in_group = var.n_variables();
2312 
2313  if (var.type().family == SCALAR &&
2314  (!elem ||
2315  var.active_on_subdomain(elem->subdomain_id())))
2316  {
2317  for (unsigned int vig=0; vig != vars_in_group; ++vig)
2318  {
2319 #ifdef DEBUG
2320  tot_size += var.type().order;
2321 #endif
2322  std::vector<dof_id_type> di_new;
2323  this->SCALAR_dof_indices(di_new,var.number(vig));
2324  di.insert( di.end(), di_new.begin(), di_new.end());
2325  }
2326  }
2327  else if (elem)
2328  for (unsigned int vig=0; vig != vars_in_group; ++vig)
2329  {
2330  _dof_indices(*elem, elem->p_level(), di, vg, vig,
2331  elem->get_nodes(), n_nodes, var.number(vig)
2332 #ifdef DEBUG
2333  , tot_size
2334 #endif
2335  );
2336  }
2337  }
2338 
2339 #ifdef DEBUG
2340  libmesh_assert_equal_to (tot_size, di.size());
2341 #endif
2342 }
2343 
2344 
2345 void DofMap::dof_indices (const Elem * const elem,
2346  std::vector<dof_id_type> & di,
2347  const unsigned int vn,
2348  int p_level) const
2349 {
2350  dof_indices(
2351  elem,
2352  di,
2353  vn,
2354  [](const Elem &,
2355  std::vector<dof_id_type> & dof_indices,
2356  const std::vector<dof_id_type> & scalar_dof_indices) {
2357  dof_indices.insert(dof_indices.end(), scalar_dof_indices.begin(), scalar_dof_indices.end());
2358  },
2359  [](const Elem &,
2360  unsigned int,
2361  unsigned int,
2362  std::vector<dof_id_type> & dof_indices,
2363  const dof_id_type dof) { dof_indices.push_back(dof); },
2364  p_level);
2365 }
2366 
2367 void DofMap::dof_indices (const Node * const node,
2368  std::vector<dof_id_type> & di) const
2369 {
2370  // We allow node==nullptr to request just SCALAR dofs
2371  // libmesh_assert(elem);
2372 
2373  // dof_indices() is a relatively light-weight function that is
2374  // called millions of times in normal codes. Therefore, it is not a
2375  // good candidate for logging, since the cost of the logging code
2376  // itself is roughly on par with the time required to call
2377  // dof_indices().
2378  // LOG_SCOPE("dof_indices(Node)", "DofMap");
2379 
2380  // Clear the DOF indices vector
2381  di.clear();
2382 
2383  const unsigned int n_var_groups = this->n_variable_groups();
2384  const unsigned int sys_num = this->sys_number();
2385 
2386  // Get the dof numbers
2387  for (unsigned int vg=0; vg<n_var_groups; vg++)
2388  {
2389  const VariableGroup & var = this->variable_group(vg);
2390  const unsigned int vars_in_group = var.n_variables();
2391 
2392  if (var.type().family == SCALAR)
2393  {
2394  for (unsigned int vig=0; vig != vars_in_group; ++vig)
2395  {
2396  std::vector<dof_id_type> di_new;
2397  this->SCALAR_dof_indices(di_new,var.number(vig));
2398  di.insert( di.end(), di_new.begin(), di_new.end());
2399  }
2400  }
2401  else
2402  {
2403  const int n_comp = node->n_comp_group(sys_num,vg);
2404  for (unsigned int vig=0; vig != vars_in_group; ++vig)
2405  {
2406  for (int i=0; i != n_comp; ++i)
2407  {
2408  const dof_id_type d =
2409  node->dof_number(sys_num, vg, vig, i, n_comp);
2410  libmesh_assert_not_equal_to
2411  (d, DofObject::invalid_id);
2412  di.push_back(d);
2413  }
2414  }
2415  }
2416  }
2417 }
2418 
2419 
2420 void DofMap::dof_indices (const Node * const node,
2421  std::vector<dof_id_type> & di,
2422  const unsigned int vn) const
2423 {
2424  if (vn == libMesh::invalid_uint)
2425  {
2426  this->dof_indices(node, di);
2427  return;
2428  }
2429 
2430  // We allow node==nullptr to request just SCALAR dofs
2431  // libmesh_assert(elem);
2432 
2433  // dof_indices() is a relatively light-weight function that is
2434  // called millions of times in normal codes. Therefore, it is not a
2435  // good candidate for logging, since the cost of the logging code
2436  // itself is roughly on par with the time required to call
2437  // dof_indices().
2438  // LOG_SCOPE("dof_indices(Node)", "DofMap");
2439 
2440  // Clear the DOF indices vector
2441  di.clear();
2442 
2443  const unsigned int sys_num = this->sys_number();
2444 
2445  // Get the dof numbers
2446  const unsigned int vg = this->_variable_group_numbers[vn];
2447  const VariableGroup & var = this->variable_group(vg);
2448 
2449  if (var.type().family == SCALAR)
2450  {
2451  std::vector<dof_id_type> di_new;
2452  this->SCALAR_dof_indices(di_new,vn);
2453  di.insert( di.end(), di_new.begin(), di_new.end());
2454  }
2455  else
2456  {
2457  const unsigned int vig = vn - var.number();
2458  const int n_comp = node->n_comp_group(sys_num,vg);
2459  for (int i=0; i != n_comp; ++i)
2460  {
2461  const dof_id_type d =
2462  node->dof_number(sys_num, vg, vig, i, n_comp);
2463  libmesh_assert_not_equal_to
2464  (d, DofObject::invalid_id);
2465  di.push_back(d);
2466  }
2467  }
2468 }
2469 
2470 
2471 void DofMap::dof_indices (const Elem & elem,
2472  unsigned int n,
2473  std::vector<dof_id_type> & di,
2474  const unsigned int vn) const
2475 {
2476  this->_node_dof_indices(elem, n, elem.node_ref(n), di, vn);
2477 }
2478 
2479 
2480 
2481 #ifdef LIBMESH_ENABLE_AMR
2482 
2483 void DofMap::old_dof_indices (const Elem & elem,
2484  unsigned int n,
2485  std::vector<dof_id_type> & di,
2486  const unsigned int vn) const
2487 {
2488  const DofObject & old_obj = elem.node_ref(n).get_old_dof_object_ref();
2489  this->_node_dof_indices(elem, n, old_obj, di, vn);
2490 }
2491 
2492 #endif // LIBMESH_ENABLE_AMR
2493 
2494 
2495 
2496 void DofMap::_node_dof_indices (const Elem & elem,
2497  unsigned int n,
2498  const DofObject & obj,
2499  std::vector<dof_id_type> & di,
2500  const unsigned int vn) const
2501 {
2502  // Half of this is a cut and paste of _dof_indices code below, but
2503  // duplication actually seems cleaner than creating a helper
2504  // function with a million arguments and hoping the compiler inlines
2505  // it properly into one of our most highly trafficked functions.
2506 
2507  LOG_SCOPE("_node_dof_indices()", "DofMap");
2508 
2509  const unsigned int sys_num = this->sys_number();
2510  const auto [vg, vig] =
2511  obj.var_to_vg_and_offset(sys_num,vn);
2512  const unsigned int n_comp = obj.n_comp_group(sys_num,vg);
2513 
2514  const VariableGroup & var = this->variable_group(vg);
2515  FEType fe_type = var.type();
2516  const bool extra_hanging_dofs =
2518 
2519  const bool add_p_level =
2520 #ifdef LIBMESH_ENABLE_AMR
2521  !_dont_p_refine.count(vg);
2522 #else
2523  false;
2524 #endif
2525 
2526  // There is a potential problem with h refinement. Imagine a
2527  // quad9 that has a linear FE on it. Then, on the hanging side,
2528  // it can falsely identify a DOF at the mid-edge node. This is why
2529  // we go through FEInterface instead of obj->n_comp() directly.
2530  const unsigned int nc =
2531  FEInterface::n_dofs_at_node(fe_type, &elem, n, add_p_level);
2532 
2533  // If this is a non-vertex on a hanging node with extra
2534  // degrees of freedom, we use the non-vertex dofs (which
2535  // come in reverse order starting from the end, to
2536  // simplify p refinement)
2537  if (extra_hanging_dofs && nc && !elem.is_vertex(n))
2538  {
2539  const int dof_offset = n_comp - nc;
2540 
2541  // We should never have fewer dofs than necessary on a
2542  // node unless we're getting indices on a parent element,
2543  // and we should never need the indices on such a node
2544  if (dof_offset < 0)
2545  {
2546  libmesh_assert(!elem.active());
2547  di.resize(di.size() + nc, DofObject::invalid_id);
2548  }
2549  else
2550  for (unsigned int i = dof_offset; i != n_comp; ++i)
2551  {
2552  const dof_id_type d =
2553  obj.dof_number(sys_num, vg, vig, i, n_comp);
2554  libmesh_assert_not_equal_to (d, DofObject::invalid_id);
2555  di.push_back(d);
2556  }
2557  }
2558  // If this is a vertex or an element without extra hanging
2559  // dofs, our dofs come in forward order coming from the
2560  // beginning. But we still might not have all those dofs, in cases
2561  // where a subdomain-restricted variable just had its subdomain
2562  // expanded.
2563  else
2564  {
2565  const unsigned int good_nc =
2566  std::min(static_cast<unsigned int>(n_comp), nc);
2567  for (unsigned int i=0; i != good_nc; ++i)
2568  {
2569  const dof_id_type d =
2570  obj.dof_number(sys_num, vg, vig, i, n_comp);
2571  libmesh_assert_not_equal_to (d, DofObject::invalid_id);
2572  di.push_back(d);
2573  }
2574  for (unsigned int i=good_nc; i != nc; ++i)
2575  di.push_back(DofObject::invalid_id);
2576  }
2577 }
2578 
2579 void
2581  int p_level,
2582  std::vector<dof_id_type> & di,
2583  const unsigned int vg,
2584  const unsigned int vig,
2585  const Node * const * nodes,
2586  unsigned int n_nodes,
2587  const unsigned int v
2588 #ifdef DEBUG
2589  ,
2590  std::size_t & tot_size
2591 #endif
2592 ) const
2593 {
2594  _dof_indices(elem,
2595  p_level,
2596  di,
2597  vg,
2598  vig,
2599  nodes,
2600  n_nodes,
2601  v,
2602 #ifdef DEBUG
2603  tot_size,
2604 #endif
2605  [](const Elem &,
2606  unsigned int,
2607  unsigned int,
2608  std::vector<dof_id_type> & functor_di,
2609  const dof_id_type dof) { functor_di.push_back(dof); });
2610 }
2611 
2612 void DofMap::SCALAR_dof_indices (std::vector<dof_id_type> & di,
2613  const unsigned int vn,
2614 #ifdef LIBMESH_ENABLE_AMR
2615  const bool old_dofs
2616 #else
2617  const bool
2618 #endif
2619  ) const
2620 {
2621  LOG_SCOPE("SCALAR_dof_indices()", "DofMap");
2622 
2623  libmesh_assert(this->variable(vn).type().family == SCALAR);
2624 
2625 #ifdef LIBMESH_ENABLE_AMR
2626  // If we're asking for old dofs then we'd better have some
2627  if (old_dofs)
2628  libmesh_assert_greater_equal(n_old_dofs(), n_SCALAR_dofs());
2629 
2630  dof_id_type my_idx = old_dofs ?
2631  this->_first_old_scalar_df[vn] : this->_first_scalar_df[vn];
2632 #else
2633  dof_id_type my_idx = this->_first_scalar_df[vn];
2634 #endif
2635 
2636  libmesh_assert_not_equal_to(my_idx, DofObject::invalid_id);
2637 
2638  // The number of SCALAR dofs comes from the variable order
2639  const int n_dofs_vn = this->variable(vn).type().order.get_order();
2640 
2641  di.resize(n_dofs_vn);
2642  for (int i = 0; i != n_dofs_vn; ++i)
2643  di[i] = my_idx++;
2644 }
2645 
2646 
2647 
2648 bool DofMap::semilocal_index (dof_id_type dof_index) const
2649 {
2650  // If it's not in the local indices
2651  if (!this->local_index(dof_index))
2652  {
2653  // and if it's not in the ghost indices, then we're not
2654  // semilocal
2655  if (!std::binary_search(_send_list.begin(), _send_list.end(), dof_index))
2656  return false;
2657  }
2658 
2659  return true;
2660 }
2661 
2662 
2663 
2664 bool DofMap::all_semilocal_indices (const std::vector<dof_id_type> & dof_indices_in) const
2665 {
2666  // We're all semilocal unless we find a counterexample
2667  for (const auto & di : dof_indices_in)
2668  if (!this->semilocal_index(di))
2669  return false;
2670 
2671  return true;
2672 }
2673 
2674 
2675 
2676 template <typename DofObjectSubclass>
2677 bool DofMap::is_evaluable(const DofObjectSubclass & obj,
2678  unsigned int var_num) const
2679 {
2680  // Everything is evaluable on a local object
2681  if (obj.processor_id() == this->processor_id())
2682  return true;
2683 
2684  std::vector<dof_id_type> di;
2685 
2686  if (var_num == libMesh::invalid_uint)
2687  this->dof_indices(&obj, di);
2688  else
2689  this->dof_indices(&obj, di, var_num);
2690 
2691  return this->all_semilocal_indices(di);
2692 }
2693 
2694 
2695 
2696 #ifdef LIBMESH_ENABLE_AMR
2697 
2698 void DofMap::old_dof_indices (const Elem * const elem,
2699  std::vector<dof_id_type> & di,
2700  const unsigned int vn) const
2701 {
2702  LOG_SCOPE("old_dof_indices()", "DofMap");
2703 
2704  libmesh_assert(elem);
2705 
2706  const ElemType type = elem->type();
2707  const unsigned int sys_num = this->sys_number();
2708  const unsigned int n_var_groups = this->n_variable_groups();
2709 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2710  const bool is_inf = elem->infinite();
2711 #endif
2712 
2713  // If we have dof indices stored on the elem, and there's no chance
2714  // that we only have those indices because we were just p refined,
2715  // then we should have old dof indices too.
2716  libmesh_assert(!elem->has_dofs(sys_num) ||
2718  elem->get_old_dof_object());
2719 
2720  // Clear the DOF indices vector.
2721  di.clear();
2722 
2723  // Determine the nodes contributing to element elem
2724  std::vector<const Node *> elem_nodes;
2725  const Node * const * nodes_ptr;
2726  unsigned int n_nodes;
2727  if (elem->type() == TRI3SUBDIVISION)
2728  {
2729  // Subdivision surface FE require the 1-ring around elem
2730  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
2731  MeshTools::Subdivision::find_one_ring(sd_elem, elem_nodes);
2732  nodes_ptr = elem_nodes.data();
2733  n_nodes = cast_int<unsigned int>(elem_nodes.size());
2734  }
2735  else
2736  {
2737  // All other FE use only the nodes of elem itself
2738  nodes_ptr = elem->get_nodes();
2739  n_nodes = elem->n_nodes();
2740  }
2741 
2742  // Get the dof numbers
2743  for (unsigned int vg=0; vg<n_var_groups; vg++)
2744  {
2745  const VariableGroup & var = this->variable_group(vg);
2746  const unsigned int vars_in_group = var.n_variables();
2747 
2748  for (unsigned int vig=0; vig<vars_in_group; vig++)
2749  {
2750  const unsigned int v = var.number(vig);
2751  if ((vn == v) || (vn == libMesh::invalid_uint))
2752  {
2753  if (var.type().family == SCALAR &&
2754  (!elem ||
2755  var.active_on_subdomain(elem->subdomain_id())))
2756  {
2757  // We asked for this variable, so add it to the vector.
2758  std::vector<dof_id_type> di_new;
2759  this->SCALAR_dof_indices(di_new,v,true);
2760  di.insert( di.end(), di_new.begin(), di_new.end());
2761  }
2762  else
2763  if (var.active_on_subdomain(elem->subdomain_id()))
2764  { // Do this for all the variables if one was not specified
2765  // or just for the specified variable
2766 
2767  FEType fe_type = var.type();
2768  const bool add_p_level =
2769 #ifdef LIBMESH_ENABLE_AMR
2770  !_dont_p_refine.count(vg);
2771 #else
2772  false;
2773 #endif
2774  // Increase the polynomial order on p refined elements,
2775  // but make sure you get the right polynomial order for
2776  // the OLD degrees of freedom
2777  int p_adjustment = 0;
2778  if (elem->p_refinement_flag() == Elem::JUST_REFINED)
2779  {
2780  libmesh_assert_greater (elem->p_level(), 0);
2781  p_adjustment = -1;
2782  }
2783  else if (elem->p_refinement_flag() == Elem::JUST_COARSENED)
2784  {
2785  p_adjustment = 1;
2786  }
2787  p_adjustment *= add_p_level;
2788 
2789  // Compute the net amount of "extra" order, including Elem::p_level()
2790  int extra_order = int(add_p_level*elem->p_level()) + p_adjustment;
2791 
2792  const bool extra_hanging_dofs =
2794 
2795  const FEInterface::n_dofs_at_node_ptr ndan =
2796  FEInterface::n_dofs_at_node_function(fe_type, elem);
2797 
2798  // Get the node-based DOF numbers
2799  for (unsigned int n=0; n<n_nodes; n++)
2800  {
2801  const Node * node = nodes_ptr[n];
2802  const DofObject & old_dof_obj = node->get_old_dof_object_ref();
2803 
2804  // There is a potential problem with h refinement. Imagine a
2805  // quad9 that has a linear FE on it. Then, on the hanging side,
2806  // it can falsely identify a DOF at the mid-edge node. This is why
2807  // we call FEInterface instead of node->n_comp() directly.
2808  const unsigned int nc =
2809 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2810  is_inf ?
2811  FEInterface::n_dofs_at_node(var.type(), extra_order, elem, n) :
2812 #endif
2813  ndan (type, var.type().order + extra_order, n);
2814 
2815  const int n_comp = old_dof_obj.n_comp_group(sys_num,vg);
2816 
2817  // If this is a non-vertex on a hanging node with extra
2818  // degrees of freedom, we use the non-vertex dofs (which
2819  // come in reverse order starting from the end, to
2820  // simplify p refinement)
2821  if (extra_hanging_dofs && !elem->is_vertex(n))
2822  {
2823  const int dof_offset = n_comp - nc;
2824 
2825  // We should never have fewer dofs than necessary on a
2826  // node unless we're getting indices on a parent element
2827  // or a just-coarsened element
2828  if (dof_offset < 0)
2829  {
2830  libmesh_assert(!elem->active() || elem->refinement_flag() ==
2832  di.resize(di.size() + nc, DofObject::invalid_id);
2833  }
2834  else
2835  for (int i=n_comp-1; i>=dof_offset; i--)
2836  {
2837  const dof_id_type d =
2838  old_dof_obj.dof_number(sys_num, vg, vig, i, n_comp);
2839 
2840  // On a newly-expanded subdomain, we
2841  // may have some DoFs that didn't
2842  // exist in the old system, in which
2843  // case we can't assert this:
2844  // libmesh_assert_not_equal_to (d, DofObject::invalid_id);
2845 
2846  di.push_back(d);
2847  }
2848  }
2849  // If this is a vertex or an element without extra hanging
2850  // dofs, our dofs come in forward order coming from the
2851  // beginning. But we still might not have all
2852  // those dofs on the old_dof_obj, in cases
2853  // where a subdomain-restricted variable just
2854  // had its subdomain expanded.
2855  else
2856  {
2857  const unsigned int old_nc =
2858  std::min(static_cast<unsigned int>(n_comp), nc);
2859  for (unsigned int i=0; i != old_nc; ++i)
2860  {
2861  const dof_id_type d =
2862  old_dof_obj.dof_number(sys_num, vg, vig, i, n_comp);
2863 
2864  libmesh_assert_not_equal_to (d, DofObject::invalid_id);
2865 
2866  di.push_back(d);
2867  }
2868  for (unsigned int i=old_nc; i != nc; ++i)
2869  di.push_back(DofObject::invalid_id);
2870  }
2871  }
2872 
2873  // If there are any element-based DOF numbers, get them
2874  const unsigned int nc =
2875  FEInterface::n_dofs_per_elem(fe_type, extra_order, elem);
2876 
2877  if (nc != 0)
2878  {
2879  const DofObject & old_dof_obj = elem->get_old_dof_object_ref();
2880 
2881  const unsigned int n_comp =
2882  old_dof_obj.n_comp_group(sys_num,vg);
2883 
2884  if (old_dof_obj.n_systems() > sys_num &&
2885  nc <= n_comp)
2886  {
2887 
2888  for (unsigned int i=0; i<nc; i++)
2889  {
2890  const dof_id_type d =
2891  old_dof_obj.dof_number(sys_num, vg, vig, i, n_comp);
2892 
2893  di.push_back(d);
2894  }
2895  }
2896  else
2897  {
2898  // We should never have fewer dofs than
2899  // necessary on an element unless we're
2900  // getting indices on a parent element, a
2901  // just-coarsened element ... or a
2902  // subdomain-restricted variable with a
2903  // just-expanded subdomain
2904  // libmesh_assert(!elem->active() || fe_type.family == LAGRANGE ||
2905  // elem->refinement_flag() == Elem::JUST_COARSENED);
2906  di.resize(di.size() + nc, DofObject::invalid_id);
2907  }
2908  }
2909  }
2910  }
2911  } // end loop over variables within group
2912  } // end loop over variable groups
2913 }
2914 
2915 #endif // LIBMESH_ENABLE_AMR
2916 
2917 
2918 #ifdef LIBMESH_ENABLE_CONSTRAINTS
2919 
2920 void DofMap::find_connected_dofs (std::vector<dof_id_type> & elem_dofs) const
2921 {
2922  typedef std::set<dof_id_type> RCSet;
2923 
2924  // First insert the DOFS we already depend on into the set.
2925  RCSet dof_set (elem_dofs.begin(), elem_dofs.end());
2926 
2927  bool done = true;
2928 
2929  // Next insert any dofs those might be constrained in terms
2930  // of. Note that in this case we may not be done: Those may
2931  // in turn depend on others. So, we need to repeat this process
2932  // in that case until the system depends only on unconstrained
2933  // degrees of freedom.
2934  for (const auto & dof : elem_dofs)
2935  if (this->is_constrained_dof(dof))
2936  {
2937  // If the DOF is constrained
2938  DofConstraints::const_iterator
2939  pos = _dof_constraints.find(dof);
2940 
2941  libmesh_assert (pos != _dof_constraints.end());
2942 
2943  const DofConstraintRow & constraint_row = pos->second;
2944 
2945  // adaptive p refinement currently gives us lots of empty constraint
2946  // rows - we should optimize those DoFs away in the future. [RHS]
2947  //libmesh_assert (!constraint_row.empty());
2948 
2949  // Add the DOFs this dof is constrained in terms of.
2950  // note that these dofs might also be constrained, so
2951  // we will need to call this function recursively.
2952  for (const auto & pr : constraint_row)
2953  if (!dof_set.count (pr.first))
2954  {
2955  dof_set.insert (pr.first);
2956  done = false;
2957  }
2958  }
2959 
2960 
2961  // If not done then we need to do more work
2962  // (obviously :-) )!
2963  if (!done)
2964  {
2965  // Fill the vector with the contents of the set
2966  elem_dofs.clear();
2967  elem_dofs.insert (elem_dofs.end(),
2968  dof_set.begin(), dof_set.end());
2969 
2970 
2971  // May need to do this recursively. It is possible
2972  // that we just replaced a constrained DOF with another
2973  // constrained DOF.
2974  this->find_connected_dofs (elem_dofs);
2975 
2976  } // end if (!done)
2977 }
2978 
2979 #endif // LIBMESH_ENABLE_CONSTRAINTS
2980 
2981 
2982 
2983 void DofMap::print_info(std::ostream & os) const
2984 {
2985  os << this->get_info();
2986 }
2987 
2988 
2989 
2990 std::string DofMap::get_info() const
2991 {
2992  std::ostringstream os;
2993 
2994  // If we didn't calculate the exact sparsity pattern, the threaded
2995  // sparsity pattern assembly may have just given us an upper bound
2996  // on sparsity.
2997  const char * may_equal = " <= ";
2998 
2999  // If we calculated the exact sparsity pattern, then we can report
3000  // exact bandwidth figures:
3001  for (const auto & mat : _matrices)
3002  if (mat->need_full_sparsity_pattern())
3003  may_equal = " = ";
3004 
3005  dof_id_type max_n_nz = 0, max_n_oz = 0;
3006  long double avg_n_nz = 0, avg_n_oz = 0;
3007 
3008  if (_sp)
3009  {
3010  for (const auto & val : _sp->get_n_nz())
3011  {
3012  max_n_nz = std::max(max_n_nz, val);
3013  avg_n_nz += val;
3014  }
3015 
3016  std::size_t n_nz_size = _sp->get_n_nz().size();
3017 
3018  this->comm().max(max_n_nz);
3019  this->comm().sum(avg_n_nz);
3020  this->comm().sum(n_nz_size);
3021 
3022  avg_n_nz /= std::max(n_nz_size,std::size_t(1));
3023 
3024  for (const auto & val : _sp->get_n_oz())
3025  {
3026  max_n_oz = std::max(max_n_oz, val);
3027  avg_n_oz += val;
3028  }
3029 
3030  std::size_t n_oz_size = _sp->get_n_oz().size();
3031 
3032  this->comm().max(max_n_oz);
3033  this->comm().sum(avg_n_oz);
3034  this->comm().sum(n_oz_size);
3035 
3036  avg_n_oz /= std::max(n_oz_size,std::size_t(1));
3037  }
3038 
3039  os << " DofMap Sparsity\n Average On-Processor Bandwidth"
3040  << may_equal << avg_n_nz << '\n';
3041 
3042  os << " Average Off-Processor Bandwidth"
3043  << may_equal << avg_n_oz << '\n';
3044 
3045  os << " Maximum On-Processor Bandwidth"
3046  << may_equal << max_n_nz << '\n';
3047 
3048  os << " Maximum Off-Processor Bandwidth"
3049  << may_equal << max_n_oz << std::endl;
3050 
3051 #ifdef LIBMESH_ENABLE_CONSTRAINTS
3052 
3053  std::size_t n_constraints = 0, max_constraint_length = 0,
3054  n_rhss = 0;
3055  long double avg_constraint_length = 0.;
3056 
3057  for (const auto & [constrained_dof, row] : _dof_constraints)
3058  {
3059  // Only count local constraints, then sum later
3060  if (!this->local_index(constrained_dof))
3061  continue;
3062 
3063  std::size_t rowsize = row.size();
3064 
3065  max_constraint_length = std::max(max_constraint_length,
3066  rowsize);
3067  avg_constraint_length += rowsize;
3068  n_constraints++;
3069 
3070  if (_primal_constraint_values.count(constrained_dof))
3071  n_rhss++;
3072  }
3073 
3074  this->comm().sum(n_constraints);
3075  this->comm().sum(n_rhss);
3076  this->comm().sum(avg_constraint_length);
3077  this->comm().max(max_constraint_length);
3078 
3079  os << " DofMap Constraints\n Number of DoF Constraints = "
3080  << n_constraints;
3081  if (n_rhss)
3082  os << '\n'
3083  << " Number of Heterogenous Constraints= " << n_rhss;
3084  if (n_constraints)
3085  {
3086  avg_constraint_length /= n_constraints;
3087 
3088  os << '\n'
3089  << " Average DoF Constraint Length= " << avg_constraint_length;
3090  }
3091 
3092 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3093  std::size_t n_node_constraints = 0, max_node_constraint_length = 0,
3094  n_node_rhss = 0;
3095  long double avg_node_constraint_length = 0.;
3096 
3097  for (const auto & [node, pr] : _node_constraints)
3098  {
3099  // Only count local constraints, then sum later
3100  if (node->processor_id() != this->processor_id())
3101  continue;
3102 
3103  const NodeConstraintRow & row = pr.first;
3104  std::size_t rowsize = row.size();
3105 
3106  max_node_constraint_length = std::max(max_node_constraint_length,
3107  rowsize);
3108  avg_node_constraint_length += rowsize;
3109  n_node_constraints++;
3110 
3111  if (pr.second != Point(0))
3112  n_node_rhss++;
3113  }
3114 
3115  this->comm().sum(n_node_constraints);
3116  this->comm().sum(n_node_rhss);
3117  this->comm().sum(avg_node_constraint_length);
3118  this->comm().max(max_node_constraint_length);
3119 
3120  os << "\n Number of Node Constraints = " << n_node_constraints;
3121  if (n_node_rhss)
3122  os << '\n'
3123  << " Number of Heterogenous Node Constraints= " << n_node_rhss;
3124  if (n_node_constraints)
3125  {
3126  avg_node_constraint_length /= n_node_constraints;
3127  os << "\n Maximum Node Constraint Length= " << max_node_constraint_length
3128  << '\n'
3129  << " Average Node Constraint Length= " << avg_node_constraint_length;
3130  }
3131 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
3132 
3133  os << std::endl;
3134 
3135 #endif // LIBMESH_ENABLE_CONSTRAINTS
3136 
3137  return os.str();
3138 }
3139 
3141 {
3142  _sc = std::make_unique<StaticCondensationDofMap>(mesh, sys, *this);
3143 }
3144 
3146 {
3147  if (_sc)
3148  _sc->reinit();
3149 }
3150 
3151 template LIBMESH_EXPORT bool DofMap::is_evaluable<Elem>(const Elem &, unsigned int) const;
3152 template LIBMESH_EXPORT bool DofMap::is_evaluable<Node>(const Node &, unsigned int) const;
3153 
3154 } // namespace libMesh
std::vector< VariableGroup > _variable_groups
The finite element type for each variable group.
Definition: dof_map.h:1957
void find_connected_dofs(std::vector< dof_id_type > &elem_dofs) const
Finds all the DOFS associated with the element DOFs elem_dofs.
Definition: dof_map.C:2920
std::unique_ptr< SparsityPattern::Build > _sp
The sparsity pattern of the global matrix.
Definition: dof_map.h:2089
static unsigned int n_dofs_per_elem(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:530
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
OStreamProxy err
virtual void clear()
Definition: dof_map_base.C:71
FEFamily family
The type of finite element.
Definition: fe_type.h:221
T command_line_next(std::string name, T default_value)
Use GetPot&#39;s search()/next() functions to get following arguments from the command line...
Definition: libmesh.C:1078
A class holding degree of freedom information pertinent to static condensation.
dof_id_type vg_dof_base(const unsigned int s, const unsigned int vg) const
VariableGroup DoF indices are indexed as id = base + var_in_vg*ncomp + comp This method allows for di...
Definition: dof_object.h:1317
bool _implicit_neighbor_dofs_initialized
Bools to indicate if we override the –implicit_neighbor_dofs commandline options.
Definition: dof_map.h:2159
ElemType
Defines an enum for geometric element types.
std::vector< GhostingFunctor * >::const_iterator GhostingFunctorIterator
Iterator type for coupling and algebraic ghosting functor ranges.
Definition: dof_map.h:315
unsigned int n_variable_groups() const
Definition: dof_map.h:632
bool _implicit_neighbor_dofs
Definition: dof_map.h:2160
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:1032
bool is_prepared() const
Definition: mesh_base.h:198
DefaultCoupling & default_coupling()
Default coupling functor.
Definition: dof_map.h:378
constraint_rows_type & get_constraint_rows()
Constraint rows accessors.
Definition: mesh_base.h:1709
A Node is like a Point, but with more information.
Definition: node.h:52
This abstract base class defines the interface by which library code and user code can report associa...
dof_id_type n_SCALAR_dofs() const
Definition: dof_map.h:685
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
Build the constraint matrix C and the forcing vector H associated with the element degree of freedom ...
~DofMap()
Destructor.
Definition: dof_map.C:202
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:310
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:1002
void reinit(MeshBase &mesh, const std::map< const Node *, std::set< subdomain_id_type >> &constraining_subdomains)
Reinitialize the underlying data structures conformal to the current mesh.
Definition: dof_map.C:503
bool _error_on_constraint_loop
This flag indicates whether or not we do an opt-mode check for the presence of constraint loops...
Definition: dof_map.h:1941
unsigned int n_var_groups(const unsigned int s) const
Definition: dof_object.h:957
void * _extra_sparsity_context
A pointer associated with the extra sparsity that can optionally be passed in.
Definition: dof_map.h:2013
void extract_local_vector(const NumericVector< Number > &Ug, const std::vector< dof_id_type > &dof_indices, DenseVectorBase< Number > &Ue) const
Builds the local element vector Ue from the global vector Ug, accounting for any constrained degrees ...
Definition: dof_map.C:2147
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2229
dof_id_type n_dofs() const
Definition: dof_map_base.h:105
The definition of the const_element_iterator struct.
Definition: mesh_base.h:2222
std::size_t distribute_dofs(MeshBase &)
Distribute dofs on the current mesh.
Definition: dof_map.C:978
void set_implicit_neighbor_dofs(bool implicit_neighbor_dofs)
Allow the implicit_neighbor_dofs flag to be set programmatically.
Definition: dof_map.C:1920
void add_default_ghosting()
Add the default functor(s) for coupling and algebraic ghosting.
Definition: dof_map.C:2024
We&#39;re using a class instead of a typedef to allow forward declarations and future flexibility...
void local_variable_indices(T &idx, const MeshBase &mesh, unsigned int var_num) const
If T == dof_id_type, counts, if T == std::vector<dof_id_type>, fills an array of, those dof indices w...
Definition: dof_map.C:1151
bool is_periodic_boundary(const boundary_id_type boundaryid) const
Definition: dof_map.C:215
std::vector< dof_id_type > _first_df
First DOF index on processor p.
Definition: dof_map_base.h:154
GhostingFunctorIterator algebraic_ghosting_functors_begin() const
Beginning of range of algebraic ghosting functors.
Definition: dof_map.h:428
virtual numeric_index_type size() const =0
virtual void set_mesh(const MeshBase *mesh)
It should be called after cloning a ghosting functor.
RefinementState p_refinement_flag() const
Definition: elem.h:3226
std::unique_ptr< StaticCondensationDofMap > _sc
Static condensation class.
Definition: dof_map.h:2175
std::vector< dof_id_type > _send_list
A list containing all the global DOF indices that affect the solution on my processor.
Definition: dof_map.h:1996
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:2230
std::unique_ptr< SparsityPattern::Build > build_sparsity(const MeshBase &mesh, bool calculate_constrained=false, bool use_condensed_system=false) const
Builds a sparsity pattern for matrices using the current degree-of-freedom numbering and coupling...
Definition: dof_map.C:61
void set_verify_dirichlet_bc_consistency(bool val)
Set the _verify_dirichlet_bc_consistency flag.
Definition: dof_map.C:1926
void sum(T &r) const
bool is_attached(SparseMatrix< Number > &matrix)
Matrices should not be attached more than once.
Definition: dof_map.C:329
void attach_matrix(SparseMatrix< Number > &matrix)
Additional matrices may be attached to this DofMap.
Definition: dof_map.C:274
void remove_ghosting_functor(GhostingFunctor &ghosting_functor)
Removes a functor which was previously added to the set of ghosting functors.
Definition: mesh_base.C:970
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
void clear_send_list()
Clears the _send_list vector.
Definition: dof_map.h:507
MeshBase & mesh
DefaultCoupling & default_algebraic_ghosting()
Default algebraic ghosting functor.
Definition: dof_map.h:440
void active_family_tree(std::vector< const Elem *> &active_family, bool reset=true) const
Same as the family_tree() member, but only adds the active children.
Definition: elem.C:2133
void attach_sparsity_pattern(const SparsityPattern::Build &sp)
Set a pointer to a sparsity pattern to use.
Definition: sparse_matrix.C:82
This proxy class acts like a container of indices from a single coupling row.
const Parallel::Communicator & comm() const
unsigned int p_level() const
Definition: elem.h:3108
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:215
unsigned int m() const
bool use_coupled_neighbor_dofs(const MeshBase &mesh) const
Tells other library functions whether or not this problem includes coupling between dofs in neighbori...
Definition: dof_map.C:1932
std::map< const Elem *, const CouplingMatrix *, CompareDofObjectsByPIDAndThenID > map_type
What elements do we care about and what variables do we care about on each element?
void reinit_static_condensation()
Calls reinit on the static condensation map if it exists.
Definition: dof_map.C:3145
The StoredRange class defines a contiguous, divisible set of objects.
Definition: stored_range.h:54
void set_vg_dof_base(const unsigned int s, const unsigned int vg, const dof_id_type db)
VariableGroup DoF indices are indexed as id = base + var_in_vg*ncomp + comp This method allows for di...
Definition: dof_object.h:1297
GhostingFunctorIterator algebraic_ghosting_functors_end() const
End of range of algebraic ghosting functors.
Definition: dof_map.h:434
virtual void augment_send_list(std::vector< dof_id_type > &send_list)=0
User-defined function to augment the send list.
The libMesh namespace provides an interface to certain functionality in the library.
bool has_dofs(const unsigned int s=libMesh::invalid_uint) const
Definition: dof_object.h:1226
virtual void zero()=0
Set every element in the vector to 0.
void set_error_on_constraint_loop(bool error_on_constraint_loop)
Definition: dof_map.C:241
Real distance(const Point &p)
void add_coupling_functor(GhostingFunctor &coupling_functor, bool to_mesh=true)
Adds a functor which can specify coupling requirements for creation of sparse matrices.
Definition: dof_map.C:2033
std::string get_info() const
Gets summary info about the sparsity bandwidth and constraints.
Definition: dof_map.C:2990
static unsigned int max_order(const FEType &fe_t, const ElemType &el_t)
unsigned int sys_number() const
Definition: dof_map.h:2182
void SCALAR_dof_indices(std::vector< dof_id_type > &di, const unsigned int vn, const bool old_dofs=false) const
Fills the vector di with the global degree of freedom indices corresponding to the SCALAR variable vn...
Definition: dof_map.C:2612
uint8_t processor_id_type
Definition: id_types.h:104
This is the MeshBase class.
Definition: mesh_base.h:75
std::vector< dof_id_type > _first_scalar_df
First DOF index for SCALAR variable v, or garbage for non-SCALAR variable v.
Definition: dof_map.h:1990
const Variable & variable(const unsigned int c) const override
Definition: dof_map.h:2200
This class implements the default algebraic coupling in libMesh: elements couple to themselves...
AdjointDofConstraintValues _adjoint_constraint_values
Definition: dof_map.h:2120
bool _constrained_sparsity_construction
This flag indicates whether or not we explicitly take constraint equations into account when computin...
Definition: dof_map.h:1947
dof_id_type end_dof() const
Definition: dof_map_base.h:83
AugmentSendList * _augment_send_list
Function object to call to add extra entries to the send list.
Definition: dof_map.h:2018
void libmesh_assert_valid_dof_ids(const MeshBase &mesh, unsigned int sysnum=libMesh::invalid_uint)
A function for verifying that degree of freedom indexing matches across processors.
Definition: mesh_tools.C:1878
bool need_full_sparsity_pattern
Default false; set to true if any attached matrix requires a full sparsity pattern.
Definition: dof_map.h:2082
void distribute_local_dofs_var_major(dof_id_type &next_free_dof, MeshBase &mesh, const std::map< const Node *, std::set< subdomain_id_type >> &constraining_subdomains)
Distributes the global degrees of freedom, for dofs on this processor.
Definition: dof_map.C:1447
virtual void clear() override
Free all new memory associated with the object, but restore its original state, with the mesh pointer...
Definition: dof_map.C:901
void _node_dof_indices(const Elem &elem, unsigned int n, const DofObject &obj, std::vector< dof_id_type > &di, const unsigned int vn) const
Helper function that implements the element-nodal versions of dof_indices and old_dof_indices.
Definition: dof_map.C:2496
void find_one_ring(const Tri3Subdivision *elem, std::vector< const Node *> &nodes)
Determines the 1-ring of element elem, and writes it to the nodes vector.
void reinit_send_list(MeshBase &mesh)
Clears the _send_list vector and then rebuilds it.
Definition: dof_map.C:1906
processor_id_type n_processors() const
const dof_id_type n_nodes
Definition: tecplot_io.C:67
dof_id_type first_dof() const
Definition: dof_map_base.h:73
std::unordered_map< unsigned int, unsigned int > _var_to_vg
A map from variable number to variable group number.
Definition: dof_map.h:1967
This class defines the notion of a variable in the system.
Definition: variable.h:50
bool has_static_condensation() const
Checks whether we have static condensation.
Definition: dof_map.h:1675
void add_neighbors_to_send_list(MeshBase &mesh)
Adds entries to the _send_list vector corresponding to DoFs on elements neighboring the current proce...
Definition: dof_map.C:1716
std::vector< dof_id_type > _first_old_scalar_df
First old DOF index for SCALAR variable v, or garbage for non-SCALAR variable v.
Definition: dof_map.h:2103
int8_t boundary_id_type
Definition: id_types.h:51
void(* _extra_sparsity_function)(SparsityPattern::Graph &, std::vector< dof_id_type > &n_nz, std::vector< dof_id_type > &n_oz, void *)
A function pointer to a function to call to add extra entries to the sparsity pattern.
Definition: dof_map.h:2006
const Node & node_ref(const unsigned int i) const
Definition: elem.h:2529
dof_id_type id() const
Definition: dof_object.h:828
dof_id_type _n_SCALAR_dofs
The total number of SCALAR dofs associated to all SCALAR variables.
Definition: dof_map.h:2095
void assert_no_nodes_missed(MeshBase &mesh)
Definition: dof_map.C:1592
void set_nonlocal_dof_objects(iterator_type objects_begin, iterator_type objects_end, MeshBase &mesh, dofobject_accessor objects)
Helper function for distributing dofs in parallel.
Definition: dof_map.C:352
unsigned int n_variables() const
Definition: variable.h:266
static bool extra_hanging_dofs(const FEType &fe_t)
virtual unsigned int n_nodes() const =0
std::set< std::unique_ptr< CouplingMatrix >, Utility::CompareUnderlying > CouplingMatricesSet
Definition: dof_map.h:1855
virtual unsigned int size() const =0
static const processor_id_type invalid_processor_id
An invalid processor_id to distinguish DoFs that have not been assigned to a processor.
Definition: dof_object.h:493
unsigned int n_vars(const unsigned int s, const unsigned int vg) const
Definition: dof_object.h:967
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:96
unsigned int n_variables() const override
Definition: dof_map.h:635
unsigned int n_systems() const
Definition: dof_object.h:937
void add_variable_group(VariableGroup var_group)
Add an unknown of order order and finite element type type to the system of equations.
Definition: dof_map.C:248
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
Definition: simple_range.h:57
This base class provides a minimal set of interfaces for satisfying user requests for...
Definition: dof_map_base.h:51
const Node *const * get_nodes() const
Definition: elem.h:2499
We&#39;re using a class instead of a typedef to allow forward declarations and future flexibility...
libmesh_assert(ctx)
bool is_constrained_dof(const dof_id_type dof) const
Definition: dof_map.h:2268
The Tri3Subdivision element is a three-noded subdivision surface shell element used in mechanics calc...
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:482
const VariableGroup & variable_group(const unsigned int c) const
Definition: dof_map.h:2190
dof_id_type n_old_dofs() const
Definition: dof_map_base.h:123
std::vector< GhostingFunctor * > _coupling_functors
The list of all GhostingFunctor objects to be used when coupling degrees of freedom in matrix sparsit...
Definition: dof_map.h:2070
DofConstraints _dof_constraints
Data structure containing DOF constraints.
Definition: dof_map.h:2116
DofMap(const unsigned int sys_number, MeshBase &mesh)
Constructor.
Definition: dof_map.C:136
CouplingMatrix * _dof_coupling
Degree of freedom coupling.
Definition: dof_map.h:1620
static void merge_ghost_functor_outputs(GhostingFunctor::map_type &elements_to_ghost, CouplingMatricesSet &temporary_coupling_matrices, const GhostingFunctorIterator &gf_begin, const GhostingFunctorIterator &gf_end, const MeshBase::const_element_iterator &elems_begin, const MeshBase::const_element_iterator &elems_end, processor_id_type p)
Definition: dof_map.C:1614
DofObject * get_old_dof_object()
Pointer accessor for previously public old_dof_object.
Definition: dof_object.h:96
bool active_on_subdomain(subdomain_id_type sid) const
Definition: variable.h:167
void create_static_condensation(MeshBase &mesh, System &system)
Add a static condensation class.
Definition: dof_map.C:3140
void print_info(std::ostream &os=libMesh::out) const
Prints summary info about the sparsity bandwidth and constraints.
Definition: dof_map.C:2983
void * _extra_send_list_context
A pointer associated with the extra send list that can optionally be passed in.
Definition: dof_map.h:2028
bool computed_sparsity_already() const
Returns true iff a sparsity pattern has already been computed.
Definition: dof_map.C:292
This class defines a logically grouped set of variables in the system.
Definition: variable.h:203
void update_sparsity_pattern(SparseMatrix< Number > &matrix) const
Additional matrices may be be temporarily initialized by this DofMap.
Definition: dof_map.C:303
void set_n_comp_group(const unsigned int s, const unsigned int vg, const unsigned int ncomp)
Sets the number of components for VariableGroup vg of system s associated with this DofObject...
Definition: dof_object.C:397
DofObject & get_old_dof_object_ref()
As above, but do not use in situations where the old_dof_object may be nullptr, since this function a...
Definition: dof_object.h:104
static n_dofs_at_node_ptr n_dofs_at_node_function(const unsigned int dim, const FEType &fe_t)
Definition: fe_interface.C:458
std::string enum_to_string(const T e)
int get_order() const
Explicitly request the order as an int.
Definition: fe_type.h:80
std::vector< SparseMatrix< Number > *> _matrices
Additional matrices handled by this object.
Definition: dof_map.h:1984
GhostingFunctorIterator coupling_functors_end() const
End of range of coupling functors.
Definition: dof_map.h:372
DofObject * elem_ptr(MeshBase &mesh, dof_id_type i) const
Definition: dof_map.C:344
std::vector< GhostingFunctor * > _algebraic_ghosting_functors
The list of all GhostingFunctor objects to be used when distributing ghosted vectors.
Definition: dof_map.h:2057
static FEContinuity get_continuity(const FEType &fe_type)
Returns the input FEType&#39;s FEContinuity based on the underlying FEFamily and potentially the Order...
void attach_dof_map(const DofMap &dof_map)
Set a pointer to the DofMap to use.
Definition: sparse_matrix.C:72
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
void distribute_scalar_dofs(dof_id_type &next_free_dof)
Definition: dof_map.C:1568
virtual numeric_index_type first_local_index() const =0
std::unique_ptr< PeriodicBoundaries > _periodic_boundaries
Data structure containing periodic boundaries.
Definition: dof_map.h:2136
virtual bool need_full_sparsity_pattern() const
void remove_coupling_functor(GhostingFunctor &coupling_functor)
Removes a functor which was previously added to the set of coupling functors, from both this DofMap a...
Definition: dof_map.C:2060
unsigned int(* n_dofs_at_node_ptr)(const ElemType, const Order, const unsigned int)
Definition: fe_interface.h:151
virtual void update_sparsity_pattern(const SparsityPattern::Graph &)
Updates the matrix sparsity pattern.
void process_constraints(MeshBase &)
Postprocesses any constrained degrees of freedom to be constrained only in terms of unconstrained dof...
subdomain_id_type subdomain_id() const
Definition: elem.h:2582
void max(const T &r, T &o, Request &req) const
void _dof_indices(const Elem &elem, int p_level, std::vector< dof_id_type > &di, const unsigned int vg, const unsigned int vig, const Node *const *nodes, unsigned int n_nodes, const unsigned int v #ifdef DEBUG, std::size_t &tot_size #endif) const
Helper function that gets the dof indices on the current element for a non-SCALAR type variable...
Definition: dof_map.C:2580
void invalidate_dofs(MeshBase &mesh) const
Invalidates all active DofObject dofs for this system.
Definition: dof_map.C:886
virtual bool is_vertex(const unsigned int i) const =0
bool _verify_dirichlet_bc_consistency
Flag which determines whether we should do some additional checking of the consistency of the Dirichl...
Definition: dof_map.h:2172
std::unique_ptr< DefaultCoupling > _default_coupling
The default coupling GhostingFunctor, used to implement standard libMesh sparsity pattern constructio...
Definition: dof_map.h:2036
OStreamProxy out
DofConstraintValueMap _primal_constraint_values
Definition: dof_map.h:2118
std::vector< Variable > _variables
The finite element type for each variable.
Definition: dof_map.h:1952
DofConstraints _stashed_dof_constraints
Definition: dof_map.h:2116
void remove_default_ghosting()
Remove any default ghosting functor(s).
Definition: dof_map.C:2016
std::vector< dof_id_type > _end_old_df
Last old DOF index (plus 1) on processor p.
Definition: dof_map_base.h:181
SparsityPattern::AugmentSparsityPattern * _augment_sparsity_pattern
Function object to call to add extra entries to the sparsity pattern.
Definition: dof_map.h:2001
std::map< dof_id_type, Real, std::less< dof_id_type >, Threads::scalable_allocator< std::pair< const dof_id_type, Real > > > DofConstraintRow
A row of the Dof constraint matrix.
Definition: dof_map.h:100
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
The DofObject defines an abstract base class for objects that have degrees of freedom associated with...
Definition: dof_object.h:54
void clear_sparsity()
Clears the sparsity pattern.
Definition: dof_map.C:2009
std::size_t compute_dof_info(dof_id_type n_local_dofs)
compute the key degree of freedom information given the local number of degrees of freedom on this pr...
dof_id_type n_local_dofs() const
Definition: dof_map_base.h:115
void parallel_reduce(const Range &range, Body &body)
Execute the provided reduction operation in parallel on the specified range.
Definition: threads_none.h:101
std::unordered_set< unsigned int > _dont_p_refine
A container of variable groups that we should not p-refine.
Definition: dof_map.h:2108
unsigned int number(unsigned int v) const
Definition: variable.h:304
std::pair< unsigned int, unsigned int > var_to_vg_and_offset(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:1200
void(* _extra_send_list_function)(std::vector< dof_id_type > &, void *)
A function pointer to a function to call to add extra entries to the send list.
Definition: dof_map.h:2023
DofObject * node_ptr(MeshBase &mesh, dof_id_type i) const
Definition: dof_map.C:337
Defines an abstract dense vector base class for use in Finite Element-type computations.
Definition: dof_map.h:73
std::map< const Node *, Real, std::less< const Node * >, Threads::scalable_allocator< std::pair< const Node *const, Real > > > NodeConstraintRow
A row of the Node constraint mapping.
Definition: dof_map.h:148
MeshBase & _mesh
The mesh that system uses.
Definition: dof_map.h:1977
void prepare_send_list()
Takes the _send_list vector (which may have duplicate entries) and sorts it.
Definition: dof_map.C:1864
bool all_semilocal_indices(const std::vector< dof_id_type > &dof_indices) const
Definition: dof_map.C:2664
bool on_command_line(std::string arg)
Definition: libmesh.C:987
virtual bool infinite() const =0
std::unique_ptr< DefaultCoupling > _default_evaluating
The default algebraic GhostingFunctor, used to implement standard libMesh send_list construction...
Definition: dof_map.h:2044
std::vector< unsigned int > _variable_group_numbers
The variable group number for each variable.
Definition: dof_map.h:1962
std::map< GhostingFunctor *, std::shared_ptr< GhostingFunctor > > _shared_functors
Hang on to references to any GhostingFunctor objects we were passed in shared_ptr form...
Definition: dof_map.h:2076
void remove_algebraic_ghosting_functor(GhostingFunctor &evaluable_functor)
Removes a functor which was previously added to the set of algebraic ghosting functors, from both this DofMap and from the underlying mesh.
Definition: dof_map.C:2117
unsigned int n() const
void compute_sparsity(const MeshBase &)
Computes the sparsity pattern for the matrices corresponding to proc_id and sends that data to Linear...
Definition: dof_map.C:1988
void add_algebraic_ghosting_functor(GhostingFunctor &evaluable_functor, bool to_mesh=true)
Adds a functor which can specify algebraic ghosting requirements for use with distributed vectors...
Definition: dof_map.C:2090
unsigned int n_comp_group(const unsigned int s, const unsigned int vg) const
Definition: dof_object.h:1015
processor_id_type processor_id() const
virtual T el(const unsigned int i) const =0
std::map< const Node *, std::set< subdomain_id_type > > calculate_constraining_subdomains()
We may have mesh constraint rows with dependent nodes in one subdomain but dependency nodes in anothe...
Definition: dof_map.C:1285
dof_id_type _n_old_dfs
Total number of degrees of freedom on old dof objects.
Definition: dof_map_base.h:171
bool active() const
Definition: elem.h:2941
void set_error_on_cyclic_constraint(bool error_on_cyclic_constraint)
Specify whether or not we perform an extra (opt-mode enabled) check for constraint loops...
Definition: dof_map.C:234
processor_id_type processor_id() const
Definition: dof_object.h:905
std::vector< dof_id_type > _first_old_df
First old DOF index on processor p.
Definition: dof_map_base.h:176
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
virtual numeric_index_type last_local_index() const =0
void ErrorVector unsigned int
Definition: adjoints_ex3.C:360
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
bool is_evaluable(const DofObjectSubclass &obj, unsigned int var_num=libMesh::invalid_uint) const
Definition: dof_map.C:2677
bool semilocal_index(dof_id_type dof_index) const
Definition: dof_map.C:2648
void distribute_local_dofs_node_major(dof_id_type &next_free_dof, MeshBase &mesh, const std::map< const Node *, std::set< subdomain_id_type >> &constraining_subdomains)
Distributes the global degrees of freedom for dofs on this processor.
Definition: dof_map.C:1319
unsigned int idx(const ElemType type, const unsigned int nx, const unsigned int i, const unsigned int j)
A useful inline function which replaces the macros used previously.
void old_dof_indices(const Elem &elem, unsigned int n, std::vector< dof_id_type > &di, const unsigned int vn) const
Appends to the vector di the old global degree of freedom indices for elem.node_ref(n), for one variable vn.
Definition: dof_map.C:2483
GhostingFunctorIterator coupling_functors_begin() const
Beginning of range of coupling functors.
Definition: dof_map.h:366
void add_ghosting_functor(GhostingFunctor &ghosting_functor)
Adds a functor which can specify ghosting requirements for use on distributed meshes.
Definition: mesh_base.C:948
uint8_t dof_id_type
Definition: id_types.h:67
std::vector< dof_id_type > _end_df
Last DOF index (plus 1) on processor p.
Definition: dof_map_base.h:159
bool local_index(dof_id_type dof_index) const
Definition: dof_map.h:846
const FEType & type() const
Definition: variable.h:144
NodeConstraints _node_constraints
Data structure containing DofObject constraints.
Definition: dof_map.h:2127
This class defines a coupling matrix.