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