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