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