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