Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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/node_range.h"
37 : #include "libmesh/numeric_vector.h"
38 : #include "libmesh/periodic_boundary_base.h"
39 : #include "libmesh/periodic_boundaries.h"
40 : #include "libmesh/sparse_matrix.h"
41 : #include "libmesh/sparsity_pattern.h"
42 : #include "libmesh/threads.h"
43 : #include "libmesh/static_condensation_dof_map.h"
44 : #include "libmesh/constraint_elimination_dof_map.h"
45 : #include "libmesh/system.h"
46 : #include "libmesh/parallel_fe_type.h"
47 :
48 : // TIMPI includes
49 : #include "timpi/parallel_implementation.h"
50 : #include "timpi/parallel_sync.h"
51 :
52 : // C++ Includes
53 : #include <algorithm> // for std::fill, std::equal_range, std::max, std::lower_bound, etc.
54 : #include <memory>
55 : #include <set>
56 : #include <sstream>
57 : #include <unordered_map>
58 :
59 : namespace libMesh
60 : {
61 :
62 : // ------------------------------------------------------------
63 : // DofMap member functions
64 : std::unique_ptr<SparsityPattern::Build>
65 4680 : DofMap::build_sparsity (const MeshBase & mesh,
66 : const bool calculate_constrained,
67 : const ReducedDofMap * const reduced) const
68 : {
69 1414 : libmesh_assert (mesh.is_prepared());
70 :
71 4242 : LOG_SCOPE("build_sparsity()", "DofMap");
72 :
73 : // Compute the sparsity structure of the global matrix. This can be
74 : // fed into a PetscMatrixBase to allocate exactly the number of nonzeros
75 : // necessary to store the matrix. This algorithm should be linear
76 : // in the (# of elements)*(# nodes per element)
77 :
78 : // We can be more efficient in the threaded sparsity pattern assembly
79 : // if we don't need the exact pattern. For some sparse matrix formats
80 : // a good upper bound will suffice.
81 :
82 : // See if we need to include sparsity pattern entries for coupling
83 : // between neighbor dofs
84 4680 : bool implicit_neighbor_dofs = this->use_coupled_neighbor_dofs(mesh);
85 :
86 : // We can compute the sparsity pattern in parallel on multiple
87 : // threads. The goal is for each thread to compute the full sparsity
88 : // pattern for a subset of elements. These sparsity patterns can
89 : // be efficiently merged in the SparsityPattern::Build::join()
90 : // method, especially if there is not too much overlap between them.
91 : // Even better, if the full sparsity pattern is not needed then
92 : // the number of nonzeros per row can be estimated from the
93 : // sparsity patterns created on each thread.
94 : auto sp = std::make_unique<SparsityPattern::Build>
95 : (*this,
96 3266 : this->_dof_coupling,
97 4680 : this->_coupling_functors,
98 : implicit_neighbor_dofs,
99 3266 : _need_full_sparsity_pattern,
100 : calculate_constrained,
101 4680 : reduced);
102 :
103 4680 : Threads::parallel_reduce (mesh.active_local_element_stored_range(), *sp);
104 :
105 4680 : sp->parallel_sync();
106 :
107 1414 : libmesh_assert_equal_to (sp->get_sparsity_pattern().size(), this->n_local_dofs());
108 :
109 : // Check to see if we have any extra stuff to add to the sparsity_pattern
110 4680 : if (_extra_sparsity_function)
111 : {
112 0 : if (_augment_sparsity_pattern)
113 : {
114 0 : libmesh_here();
115 0 : libMesh::out << "WARNING: You have specified both an extra sparsity function and object.\n"
116 0 : << " Are you sure this is what you meant to do??"
117 0 : << std::endl;
118 : }
119 :
120 0 : sp->apply_extra_sparsity_function(_extra_sparsity_function,
121 0 : _extra_sparsity_context);
122 : }
123 :
124 4680 : if (_augment_sparsity_pattern)
125 0 : sp->apply_extra_sparsity_object(*_augment_sparsity_pattern);
126 :
127 6094 : return sp;
128 0 : }
129 :
130 :
131 :
132 35508 : DofMap::DofMap(const unsigned int number,
133 35508 : MeshBase & mesh) :
134 : DofMapBase (mesh.comm()),
135 15128 : _dof_coupling(nullptr),
136 15128 : _error_on_constraint_loop(false),
137 15128 : _constrained_sparsity_construction(false),
138 15128 : _variables(),
139 15128 : _variable_groups(),
140 : _variable_group_numbers(),
141 15128 : _sys_number(number),
142 15128 : _mesh(mesh),
143 : _matrices(),
144 : _first_scalar_df(),
145 : _send_list(),
146 15128 : _augment_sparsity_pattern(nullptr),
147 15128 : _extra_sparsity_function(nullptr),
148 15128 : _extra_sparsity_context(nullptr),
149 15128 : _augment_send_list(nullptr),
150 15128 : _extra_send_list_function(nullptr),
151 15128 : _extra_send_list_context(nullptr),
152 15128 : _default_coupling(std::make_unique<DefaultCoupling>()),
153 15128 : _default_evaluating(std::make_unique<DefaultCoupling>()),
154 15128 : _need_full_sparsity_pattern(false),
155 15128 : _need_ghost_constraints(false),
156 15128 : _n_SCALAR_dofs(0)
157 : #ifdef LIBMESH_ENABLE_AMR
158 : , _first_old_scalar_df()
159 : #endif
160 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
161 : , _dof_constraints()
162 : , _stashed_dof_constraints()
163 : , _primal_constraint_values()
164 : , _adjoint_constraint_values()
165 : #endif
166 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
167 : , _node_constraints()
168 : #endif
169 : #ifdef LIBMESH_ENABLE_PERIODIC
170 15128 : , _periodic_boundaries(std::make_unique<PeriodicBoundaries>())
171 : #endif
172 : #ifdef LIBMESH_ENABLE_DIRICHLET
173 15128 : , _dirichlet_boundaries(std::make_unique<DirichletBoundaries>())
174 15128 : , _adjoint_dirichlet_boundaries()
175 : #endif
176 15128 : , _implicit_neighbor_dofs_initialized(false),
177 15128 : _implicit_neighbor_dofs(false),
178 15128 : _verify_dirichlet_bc_consistency(true),
179 142032 : _sc(nullptr)
180 : {
181 10190 : _matrices.clear();
182 :
183 35508 : _default_coupling->set_mesh(&_mesh);
184 35508 : _default_evaluating->set_mesh(&_mesh);
185 10190 : _default_evaluating->set_n_levels(1);
186 :
187 : #ifdef LIBMESH_ENABLE_PERIODIC
188 45698 : _default_coupling->set_periodic_boundaries(_periodic_boundaries.get());
189 45698 : _default_evaluating->set_periodic_boundaries(_periodic_boundaries.get());
190 : #endif
191 :
192 35508 : this->add_coupling_functor(*_default_coupling);
193 35508 : this->add_algebraic_ghosting_functor(*_default_evaluating);
194 35508 : }
195 :
196 :
197 :
198 : // Destructor
199 91396 : DofMap::~DofMap()
200 : {
201 35508 : this->clear();
202 :
203 : // clear() resets all but the default DofMap-based functors. We
204 : // need to remove those from the mesh too before we die.
205 45698 : _mesh.remove_ghosting_functor(*_default_coupling);
206 45698 : _mesh.remove_ghosting_functor(*_default_evaluating);
207 131528 : }
208 :
209 :
210 : #ifdef LIBMESH_ENABLE_PERIODIC
211 :
212 0 : bool DofMap::is_periodic_boundary (const boundary_id_type boundaryid) const
213 : {
214 0 : if (_periodic_boundaries->count(boundaryid) != 0)
215 0 : return true;
216 :
217 0 : return false;
218 : }
219 :
220 : #endif
221 :
222 0 : void DofMap::set_error_on_cyclic_constraint(bool error_on_cyclic_constraint)
223 : {
224 : // This function will eventually be officially libmesh_deprecated();
225 : // Call DofMap::set_error_on_constraint_loop() instead.
226 0 : set_error_on_constraint_loop(error_on_cyclic_constraint);
227 0 : }
228 :
229 7 : void DofMap::set_error_on_constraint_loop(bool error_on_constraint_loop)
230 : {
231 7 : _error_on_constraint_loop = error_on_constraint_loop;
232 7 : }
233 :
234 :
235 2285 : void DofMap::attach_matrix (SparseMatrix<Number> & matrix)
236 : {
237 674 : parallel_object_only();
238 :
239 : // We shouldn't be trying to re-attach the same matrices repeatedly
240 674 : libmesh_assert (std::find(_matrices.begin(), _matrices.end(),
241 : &matrix) == _matrices.end());
242 :
243 2285 : _matrices.push_back(&matrix);
244 :
245 2285 : this->update_sparsity_pattern(matrix);
246 :
247 2285 : if (matrix.need_full_sparsity_pattern())
248 0 : _need_full_sparsity_pattern = true;
249 2285 : }
250 :
251 :
252 :
253 2329 : bool DofMap::computed_sparsity_already() const
254 : {
255 2381 : bool computed_sparsity_already = _sp &&
256 52 : (!_sp->get_n_nz().empty() ||
257 2329 : !_sp->get_n_oz().empty());
258 2329 : this->comm().max(computed_sparsity_already);
259 2329 : return computed_sparsity_already;
260 : }
261 :
262 :
263 :
264 2285 : void DofMap::update_sparsity_pattern(SparseMatrix<Number> & matrix) const
265 : {
266 2285 : matrix.attach_dof_map (*this);
267 :
268 : // If we've already computed sparsity, then it's too late
269 : // to wait for "compute_sparsity" to help with sparse matrix
270 : // initialization, and we need to handle this matrix individually
271 2285 : if (this->computed_sparsity_already())
272 : {
273 52 : libmesh_assert(_sp.get());
274 :
275 176 : if (matrix.need_full_sparsity_pattern())
276 : {
277 : // We'd better have already computed the full sparsity
278 : // pattern if we need it here
279 0 : libmesh_assert(_need_full_sparsity_pattern);
280 :
281 0 : matrix.update_sparsity_pattern (_sp->get_sparsity_pattern());
282 : }
283 :
284 176 : matrix.attach_sparsity_pattern(*_sp);
285 : }
286 2285 : }
287 :
288 :
289 :
290 2115 : bool DofMap::is_attached (SparseMatrix<Number> & matrix)
291 : {
292 2115 : return (std::find(_matrices.begin(), _matrices.end(),
293 2737 : &matrix) != _matrices.end());
294 : }
295 :
296 :
297 :
298 12172282 : DofObject * DofMap::node_ptr(MeshBase & mesh, dof_id_type i) const
299 : {
300 12172282 : return mesh.node_ptr(i);
301 : }
302 :
303 :
304 :
305 6190956 : DofObject * DofMap::elem_ptr(MeshBase & mesh, dof_id_type i) const
306 : {
307 6190956 : return mesh.elem_ptr(i);
308 : }
309 :
310 :
311 :
312 : template <typename iterator_type>
313 67072 : void DofMap::set_nonlocal_dof_objects(iterator_type objects_begin,
314 : iterator_type objects_end,
315 : MeshBase & mesh,
316 : dofobject_accessor objects)
317 : {
318 : // This function must be run on all processors at once
319 22432 : parallel_object_only();
320 :
321 : // First, iterate over local objects to find out how many
322 : // are on each processor
323 44864 : std::unordered_map<processor_id_type, dof_id_type> ghost_objects_from_proc;
324 :
325 44864 : iterator_type it = objects_begin;
326 :
327 12456732 : for (; it != objects_end; ++it)
328 : {
329 9181619 : DofObject * obj = *it;
330 :
331 2986789 : if (obj)
332 : {
333 9181619 : processor_id_type obj_procid = obj->processor_id();
334 : // We'd better be completely partitioned by now
335 2986789 : libmesh_assert_not_equal_to (obj_procid, DofObject::invalid_processor_id);
336 9181619 : ghost_objects_from_proc[obj_procid]++;
337 : }
338 : }
339 :
340 : // Request sets to send to each processor
341 : std::map<processor_id_type, std::vector<dof_id_type>>
342 44864 : requested_ids;
343 :
344 : // We know how many of our objects live on each processor, so
345 : // reserve() space for requests from each.
346 174354 : for (auto [p, size] : ghost_objects_from_proc)
347 : {
348 143192 : if (p != this->processor_id())
349 53635 : requested_ids[p].reserve(size);
350 : }
351 :
352 12456732 : for (it = objects_begin; it != objects_end; ++it)
353 : {
354 9181619 : DofObject * obj = *it;
355 9181619 : if (obj->processor_id() != DofObject::invalid_processor_id)
356 9181619 : requested_ids[obj->processor_id()].push_back(obj->id());
357 : }
358 : #ifdef DEBUG
359 67296 : for (auto p : make_range(this->n_processors()))
360 : {
361 44864 : if (ghost_objects_from_proc.count(p))
362 35910 : libmesh_assert_equal_to (requested_ids[p].size(), ghost_objects_from_proc[p]);
363 : else
364 8954 : libmesh_assert(!requested_ids.count(p));
365 : }
366 : #endif
367 :
368 : typedef std::vector<dof_id_type> datum;
369 :
370 9135169 : auto gather_functor =
371 35462 : [this, &mesh, &objects]
372 : (processor_id_type,
373 : const std::vector<dof_id_type> & ids,
374 71820 : std::vector<datum> & data)
375 : {
376 : // Fill those requests
377 : const unsigned int
378 71820 : sys_num = this->sys_number(),
379 35910 : n_var_groups = this->n_variable_groups();
380 :
381 71820 : const std::size_t query_size = ids.size();
382 :
383 107282 : data.resize(query_size);
384 9288901 : for (auto & d : data)
385 9181619 : d.resize(2 * n_var_groups);
386 :
387 9288901 : for (std::size_t i=0; i != query_size; ++i)
388 : {
389 9181619 : DofObject * requested = (this->*objects)(mesh, ids[i]);
390 2986789 : libmesh_assert(requested);
391 2986789 : libmesh_assert_equal_to (requested->processor_id(), this->processor_id());
392 2986789 : libmesh_assert_equal_to (requested->n_var_groups(sys_num), n_var_groups);
393 19968392 : for (unsigned int vg=0; vg != n_var_groups; ++vg)
394 : {
395 : unsigned int n_comp_g =
396 3521847 : requested->n_comp_group(sys_num, vg);
397 14308620 : data[i][vg] = n_comp_g;
398 10786773 : dof_id_type my_first_dof = n_comp_g ?
399 1402019 : requested->vg_dof_base(sys_num, vg) : 0;
400 3521847 : libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
401 14308620 : data[i][n_var_groups+vg] = my_first_dof;
402 : }
403 : }
404 : };
405 :
406 9135169 : auto action_functor =
407 35462 : [this, &mesh, &objects]
408 : (processor_id_type libmesh_dbg_var(pid),
409 : const std::vector<dof_id_type> & ids,
410 71820 : const std::vector<datum> & data)
411 : {
412 : const unsigned int
413 71820 : sys_num = this->sys_number(),
414 35910 : n_var_groups = this->n_variable_groups();
415 :
416 : // Copy the id changes we've now been informed of
417 9288901 : for (auto i : index_range(ids))
418 : {
419 9181619 : DofObject * requested = (this->*objects)(mesh, ids[i]);
420 2986789 : libmesh_assert(requested);
421 2986789 : libmesh_assert_equal_to (requested->processor_id(), pid);
422 19968392 : for (unsigned int vg=0; vg != n_var_groups; ++vg)
423 : {
424 : unsigned int n_comp_g =
425 14308620 : cast_int<unsigned int>(data[i][vg]);
426 10786773 : requested->set_n_comp_group(sys_num, vg, n_comp_g);
427 10786773 : if (n_comp_g)
428 : {
429 7018077 : dof_id_type my_first_dof = data[i][n_var_groups+vg];
430 1402019 : libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
431 : requested->set_vg_dof_base
432 1402019 : (sys_num, vg, my_first_dof);
433 : }
434 : }
435 : }
436 : };
437 :
438 22432 : datum * ex = nullptr;
439 : Parallel::pull_parallel_vector_data
440 67072 : (this->comm(), requested_ids, gather_functor, action_functor, ex);
441 :
442 : #ifdef DEBUG
443 : // Double check for invalid dofs
444 3009221 : for (it = objects_begin; it != objects_end; ++it)
445 : {
446 2986789 : DofObject * obj = *it;
447 2986789 : libmesh_assert (obj);
448 2986789 : unsigned int num_variables = obj->n_vars(this->sys_number());
449 9931278 : for (unsigned int v=0; v != num_variables; ++v)
450 : {
451 : unsigned int n_comp =
452 6944489 : obj->n_comp(this->sys_number(), v);
453 6944489 : dof_id_type my_first_dof = n_comp ?
454 2745985 : obj->dof_number(this->sys_number(), v, 0) : 0;
455 6944489 : libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
456 : }
457 : }
458 : #endif
459 67072 : }
460 :
461 :
462 :
463 38844 : void DofMap::reinit
464 : (MeshBase & mesh,
465 : const std::map<const Node *, std::set<subdomain_id_type>> &
466 : constraining_subdomains)
467 : {
468 11218 : libmesh_assert (mesh.is_prepared());
469 :
470 33654 : LOG_SCOPE("reinit()", "DofMap");
471 :
472 : // This is the common case and we want to optimize for it
473 : const bool constraining_subdomains_empty =
474 11218 : constraining_subdomains.empty();
475 :
476 : // We ought to reconfigure our default coupling functor.
477 : //
478 : // The user might have removed it from our coupling functors set,
479 : // but if so, who cares, this reconfiguration is cheap.
480 :
481 : // Avoid calling set_dof_coupling() with an empty/non-nullptr
482 : // _dof_coupling matrix which may happen when there are actually no
483 : // variables on the system.
484 38844 : if (this->_dof_coupling && this->_dof_coupling->empty() && !this->n_variables())
485 0 : this->_dof_coupling = nullptr;
486 38844 : _default_coupling->set_dof_coupling(this->_dof_coupling);
487 :
488 : // By default we may want 0 or 1 levels of coupling
489 : unsigned int standard_n_levels =
490 38844 : this->use_coupled_neighbor_dofs(mesh);
491 : _default_coupling->set_n_levels
492 53001 : (std::max(_default_coupling->n_levels(), standard_n_levels));
493 :
494 : // But we *don't* want to restrict to a CouplingMatrix unless the
495 : // user does so manually; the original libMesh behavior was to put
496 : // ghost indices on the send_list regardless of variable.
497 : //_default_evaluating->set_dof_coupling(this->_dof_coupling);
498 :
499 : const unsigned int
500 22436 : sys_num = this->sys_number(),
501 11218 : n_var_groups = this->n_variable_groups();
502 :
503 : // The DofObjects need to know how many variable groups we have, and
504 : // how many variables there are in each group.
505 50062 : std::vector<unsigned int> n_vars_per_group; /**/ n_vars_per_group.reserve (n_var_groups);
506 :
507 78769 : for (unsigned int vg=0; vg<n_var_groups; vg++)
508 39925 : n_vars_per_group.push_back (this->variable_group(vg).n_variables());
509 :
510 : #ifdef LIBMESH_ENABLE_AMR
511 :
512 : //------------------------------------------------------------
513 : // Clear the old_dof_objects for all the nodes
514 : // and elements so that we can overwrite them
515 : Threads::parallel_for_each
516 38844 : (mesh.node_stored_range(),
517 1939221 : [](Node * node)
518 : {
519 6190535 : node->clear_old_dof_object();
520 1939221 : libmesh_assert (!node->get_old_dof_object());
521 1939221 : });
522 :
523 : Threads::parallel_for_each
524 38844 : (mesh.element_stored_range(),
525 1047650 : [](Elem * elem)
526 : {
527 3025363 : elem->clear_old_dof_object();
528 1047650 : libmesh_assert (!elem->get_old_dof_object());
529 1047650 : });
530 :
531 : //------------------------------------------------------------
532 : // Set the old_dof_objects for the elements that
533 : // weren't just created, if these old dof objects
534 : // had variables
535 4957312 : for (auto & elem : mesh.element_ptr_range())
536 : {
537 : // Skip the elements that were just refined
538 3493071 : if (elem->refinement_flag() == Elem::JUST_REFINED)
539 325438 : continue;
540 :
541 23654336 : for (Node & node : elem->node_ref_range())
542 20631727 : if (node.get_old_dof_object() == nullptr)
543 15697121 : if (node.has_dofs(sys_num))
544 1775611 : node.set_old_dof_object();
545 :
546 902626 : libmesh_assert (!elem->get_old_dof_object());
547 :
548 3022609 : if (elem->has_dofs(sys_num))
549 1317684 : elem->set_old_dof_object();
550 16408 : }
551 :
552 : #endif // #ifdef LIBMESH_ENABLE_AMR
553 :
554 :
555 : //------------------------------------------------------------
556 : // Then set the number of variables for each \p DofObject
557 : // equal to n_variables() for this system. This will
558 : // handle new \p DofObjects that may have just been created
559 :
560 : // All the nodes
561 : Threads::parallel_for_each
562 50062 : (mesh.node_stored_range(),
563 3922476 : [sys_num, n_vars_per_group](Node * node)
564 : {
565 6848241 : node->set_n_vars_per_group(sys_num, n_vars_per_group);
566 1939221 : });
567 :
568 : // All the elements
569 : Threads::parallel_for_each
570 38846 : (mesh.element_stored_range(),
571 2122929 : [sys_num, n_vars_per_group](Elem * elem)
572 : {
573 3493071 : elem->set_n_vars_per_group(sys_num, n_vars_per_group);
574 1047650 : });
575 :
576 : // Zero _n_SCALAR_dofs, it will be updated below.
577 38844 : this->_n_SCALAR_dofs = 0;
578 :
579 : //------------------------------------------------------------
580 : // Next allocate space for the DOF indices
581 78762 : for (unsigned int vg=0; vg<n_var_groups; vg++)
582 : {
583 11526 : const VariableGroup & vg_description = this->variable_group(vg);
584 :
585 11526 : const unsigned int n_var_in_group = vg_description.n_variables();
586 11526 : const FEType & base_fe_type = vg_description.type();
587 :
588 39925 : const bool add_p_level = base_fe_type.p_refinement;
589 :
590 : // Don't need to loop over elements for a SCALAR variable
591 : // Just increment _n_SCALAR_dofs
592 39925 : if (base_fe_type.family == SCALAR)
593 : {
594 141 : this->_n_SCALAR_dofs += base_fe_type.order.get_order()*n_var_in_group;
595 141 : continue;
596 : }
597 :
598 : // This should be constant even on p-refined elements
599 : const bool extra_hanging_dofs =
600 39784 : FEInterface::extra_hanging_dofs(base_fe_type);
601 :
602 : // For all the active elements, count vertex degrees of freedom.
603 4292066 : for (auto & elem : mesh.active_element_ptr_range())
604 : {
605 881430 : libmesh_assert(elem);
606 :
607 : // Only number dofs connected to active elements on this
608 : // processor and only for variables which are active on on
609 : // this element's subdomain or which are active on the
610 : // subdomain of a node constrained by this node.
611 : const bool active_on_elem =
612 2993423 : vg_description.active_on_subdomain(elem->subdomain_id());
613 :
614 : // If there's no way we're active on this element then we're
615 : // done
616 2993423 : if (!active_on_elem && constraining_subdomains_empty)
617 7068 : continue;
618 :
619 2986355 : FEType fe_type = base_fe_type;
620 :
621 2986355 : const ElemType type = elem->type();
622 :
623 2986383 : libmesh_error_msg_if(base_fe_type.order.get_order() >
624 : int(FEInterface::max_order(base_fe_type,type)),
625 : "ERROR: Finite element "
626 : << Utility::enum_to_string(base_fe_type.family)
627 : << " on geometric element "
628 : << Utility::enum_to_string(type)
629 : << "\nonly supports FEInterface::max_order = "
630 : << FEInterface::max_order(base_fe_type,type)
631 : << ", not fe_type.order = "
632 : << base_fe_type.order);
633 :
634 : #ifdef LIBMESH_ENABLE_AMR
635 : // Make sure we haven't done more p refinement than we can
636 : // handle
637 2986348 : if (base_fe_type.order + add_p_level*elem->p_level() >
638 : FEInterface::max_order(base_fe_type, type))
639 : {
640 : # ifdef DEBUG
641 0 : libMesh::err << "WARNING: Finite element "
642 0 : << Utility::enum_to_string(base_fe_type.family)
643 0 : << " on geometric element "
644 0 : << Utility::enum_to_string(type) << std::endl
645 0 : << "could not be p refined past FEInterface::max_order = "
646 0 : << FEInterface::max_order(base_fe_type,type)
647 0 : << std::endl;
648 : # endif
649 0 : elem->set_p_level(int(FEInterface::max_order(base_fe_type,type))
650 0 : - int(base_fe_type.order));
651 : }
652 : #endif
653 :
654 : // Allocate the vertex DOFs
655 27361350 : for (auto n : elem->node_index_range())
656 : {
657 23495614 : Node & node = elem->node_ref(n);
658 :
659 : // If we're active on the element then we're active on
660 : // its nodes. If we're not then we might *still* be
661 : // active on particular constraining nodes.
662 6597784 : bool active_on_node = active_on_elem;
663 23495614 : if (!active_on_node)
664 0 : if (auto it = constraining_subdomains.find(&node);
665 0 : it != constraining_subdomains.end())
666 0 : for (auto s : it->second)
667 0 : if (vg_description.active_on_subdomain(s))
668 : {
669 0 : active_on_node = true;
670 0 : break;
671 : }
672 :
673 13195568 : if (!active_on_node)
674 0 : continue;
675 :
676 23495614 : if (elem->is_vertex(n))
677 : {
678 : const unsigned int old_node_dofs =
679 12850236 : node.n_comp_group(sys_num, vg);
680 :
681 : const unsigned int vertex_dofs =
682 12850239 : std::max(FEInterface::n_dofs_at_node(fe_type, elem, n, add_p_level),
683 12850236 : old_node_dofs);
684 :
685 : // Some discontinuous FEs have no vertex dofs
686 12850236 : if (vertex_dofs > old_node_dofs)
687 : {
688 1898055 : node.set_n_comp_group(sys_num, vg,
689 : vertex_dofs);
690 :
691 : // Abusing dof_number to set a "this is a
692 : // vertex" flag
693 1094286 : node.set_vg_dof_base(sys_num, vg,
694 : vertex_dofs);
695 :
696 : // libMesh::out << "sys_num,vg,old_node_dofs,vertex_dofs="
697 : // << sys_num << ","
698 : // << vg << ","
699 : // << old_node_dofs << ","
700 : // << vertex_dofs << '\n',
701 : // node.debug_buffer();
702 :
703 : // libmesh_assert_equal_to (vertex_dofs, node.n_comp(sys_num, vg));
704 : // libmesh_assert_equal_to (vertex_dofs, node.vg_dof_base(sys_num, vg));
705 : }
706 : }
707 : }
708 16816 : } // done counting vertex dofs
709 :
710 : // count edge & face dofs next
711 4292048 : for (auto & elem : mesh.active_element_ptr_range())
712 : {
713 881428 : libmesh_assert(elem);
714 :
715 : // Only number dofs connected to active elements on this
716 : // processor and only for variables which are active on on
717 : // this element's subdomain or which are active on the
718 : // subdomain of a node constrained by this node.
719 : const bool active_on_elem =
720 2993416 : vg_description.active_on_subdomain(elem->subdomain_id());
721 :
722 : // If there's no way we're active on this element then we're
723 : // done
724 2993416 : if (!active_on_elem && constraining_subdomains_empty)
725 5028 : continue;
726 :
727 : // Allocate the edge and face DOFs
728 26481962 : for (auto n : elem->node_index_range())
729 : {
730 23495614 : Node & node = elem->node_ref(n);
731 :
732 : // If we're active on the element then we're active on
733 : // its nodes. If we're not then we might *still* be
734 : // active on particular constraining nodes.
735 6597784 : bool active_on_node = active_on_elem;
736 23495614 : if (!active_on_node)
737 0 : if (auto it = constraining_subdomains.find(&node);
738 0 : it != constraining_subdomains.end())
739 0 : for (auto s : it->second)
740 0 : if (vg_description.active_on_subdomain(s))
741 : {
742 0 : active_on_node = true;
743 0 : break;
744 : }
745 :
746 13195568 : if (!active_on_node)
747 0 : continue;
748 :
749 : const unsigned int old_node_dofs =
750 6597784 : node.n_comp_group(sys_num, vg);
751 :
752 20494770 : const unsigned int vertex_dofs = old_node_dofs?
753 6597784 : cast_int<unsigned int>(node.vg_dof_base (sys_num,vg)):0;
754 :
755 : const unsigned int new_node_dofs =
756 23495614 : FEInterface::n_dofs_at_node(base_fe_type, elem, n, add_p_level);
757 :
758 : // We've already allocated vertex DOFs
759 23495614 : if (elem->is_vertex(n))
760 : {
761 3764322 : libmesh_assert_greater_equal (old_node_dofs, vertex_dofs);
762 : // //if (vertex_dofs < new_node_dofs)
763 : // libMesh::out << "sys_num,vg,old_node_dofs,vertex_dofs,new_node_dofs="
764 : // << sys_num << ","
765 : // << vg << ","
766 : // << old_node_dofs << ","
767 : // << vertex_dofs << ","
768 : // << new_node_dofs << '\n',
769 : // node.debug_buffer();
770 :
771 3764322 : libmesh_assert_greater_equal (vertex_dofs, new_node_dofs);
772 : }
773 : // We need to allocate the rest
774 : else
775 : {
776 : // If this has no dofs yet, it needs no vertex
777 : // dofs, so we just give it edge or face dofs
778 10645378 : if (!old_node_dofs)
779 : {
780 8891157 : node.set_n_comp_group(sys_num, vg,
781 : new_node_dofs);
782 : // Abusing dof_number to set a "this has no
783 : // vertex dofs" flag
784 8891157 : if (new_node_dofs)
785 583470 : node.set_vg_dof_base(sys_num, vg, 0);
786 : }
787 :
788 : // If this has dofs, but has no vertex dofs,
789 : // it may still need more edge or face dofs if
790 : // we're p-refined.
791 1754221 : else if (vertex_dofs == 0)
792 : {
793 1707089 : if (new_node_dofs > old_node_dofs)
794 : {
795 63 : node.set_n_comp_group(sys_num, vg,
796 : new_node_dofs);
797 :
798 18 : node.set_vg_dof_base(sys_num, vg,
799 : vertex_dofs);
800 : }
801 : }
802 : // If this is another element's vertex,
803 : // add more (non-overlapping) edge/face dofs if
804 : // necessary
805 47132 : else if (extra_hanging_dofs)
806 : {
807 31647 : if (new_node_dofs > old_node_dofs - vertex_dofs)
808 : {
809 28933 : node.set_n_comp_group(sys_num, vg,
810 : vertex_dofs + new_node_dofs);
811 :
812 18832 : node.set_vg_dof_base(sys_num, vg,
813 : vertex_dofs);
814 : }
815 : }
816 : // If this is another element's vertex, add any
817 : // (overlapping) edge/face dofs if necessary
818 : else
819 : {
820 5076 : libmesh_assert_greater_equal (old_node_dofs, vertex_dofs);
821 15485 : if (new_node_dofs > old_node_dofs)
822 : {
823 0 : node.set_n_comp_group(sys_num, vg,
824 : new_node_dofs);
825 :
826 0 : node.set_vg_dof_base (sys_num, vg,
827 : vertex_dofs);
828 : }
829 : }
830 : }
831 : }
832 : // Allocate the element DOFs
833 : const unsigned int dofs_per_elem =
834 2986348 : FEInterface::n_dofs_per_elem(base_fe_type, elem, add_p_level);
835 :
836 2986348 : elem->set_n_comp_group(sys_num, vg, dofs_per_elem);
837 :
838 16813 : }
839 : } // end loop over variable groups
840 :
841 : // Calling DofMap::reinit() by itself makes little sense,
842 : // so we won't bother with nonlocal DofObjects.
843 : // Those will be fixed by distribute_dofs
844 :
845 : //------------------------------------------------------------
846 : // Finally, clear all the current DOF indices
847 : // (distribute_dofs expects them cleared!)
848 38837 : this->invalidate_dofs(mesh);
849 71653 : }
850 :
851 :
852 :
853 77674 : void DofMap::invalidate_dofs(MeshBase & mesh) const
854 : {
855 44864 : const unsigned int sys_num = this->sys_number();
856 :
857 : // All the nodes
858 : Threads::parallel_for_each
859 77674 : (mesh.node_stored_range(),
860 13696132 : [sys_num](Node * node)
861 : {
862 13696132 : node->invalidate_dofs(sys_num);
863 3878342 : });
864 :
865 : // All the active elements.
866 : Threads::parallel_for_each
867 77674 : (mesh.element_stored_range(),
868 6515374 : [sys_num](Elem * elem)
869 : {
870 2095236 : if (elem->active())
871 5478160 : elem->invalidate_dofs(sys_num);
872 4890682 : });
873 77674 : }
874 :
875 :
876 :
877 35602 : void DofMap::clear()
878 : {
879 35602 : DofMapBase::clear();
880 :
881 : // we don't want to clear
882 : // the coupling matrix!
883 : // It should not change...
884 : //_dof_coupling->clear();
885 : //
886 : // But it would be inconsistent to leave our coupling settings
887 : // through a clear()...
888 35602 : _dof_coupling = nullptr;
889 :
890 : // Reset ghosting functor statuses
891 : {
892 71215 : for (const auto & gf : _coupling_functors)
893 : {
894 10224 : libmesh_assert(gf);
895 35613 : _mesh.remove_ghosting_functor(*gf);
896 : }
897 10220 : this->_coupling_functors.clear();
898 :
899 : // Go back to default coupling
900 :
901 35602 : _default_coupling->set_dof_coupling(this->_dof_coupling);
902 35602 : _default_coupling->set_n_levels(this->use_coupled_neighbor_dofs(this->_mesh));
903 :
904 35602 : this->add_coupling_functor(*_default_coupling);
905 : }
906 :
907 :
908 : {
909 71243 : for (const auto & gf : _algebraic_ghosting_functors)
910 : {
911 10232 : libmesh_assert(gf);
912 35641 : _mesh.remove_ghosting_functor(*gf);
913 : }
914 10220 : this->_algebraic_ghosting_functors.clear();
915 :
916 : // Go back to default send_list generation
917 :
918 : // _default_evaluating->set_dof_coupling(this->_dof_coupling);
919 10220 : _default_evaluating->set_n_levels(1);
920 35602 : this->add_algebraic_ghosting_functor(*_default_evaluating);
921 : }
922 :
923 10220 : this->_shared_functors.clear();
924 :
925 25382 : _variables.clear();
926 25382 : _variable_groups.clear();
927 10220 : _var_to_vg.clear();
928 10220 : _variable_group_numbers.clear();
929 10220 : _array_variables.clear();
930 10220 : _first_scalar_df.clear();
931 10220 : this->clear_send_list();
932 35602 : this->clear_sparsity();
933 35602 : _need_full_sparsity_pattern = false;
934 35602 : _need_ghost_constraints = false;
935 :
936 : #ifdef LIBMESH_ENABLE_AMR
937 :
938 10220 : _dof_constraints.clear();
939 10220 : _stashed_dof_constraints.clear();
940 10220 : _primal_constraint_values.clear();
941 10220 : _adjoint_constraint_values.clear();
942 35602 : _n_old_dfs = 0;
943 10220 : _first_old_df.clear();
944 10220 : _end_old_df.clear();
945 10220 : _first_old_scalar_df.clear();
946 :
947 : #endif
948 :
949 10220 : _matrices.clear();
950 35602 : if (_sc)
951 54 : _sc->clear();
952 35602 : if (_constraint_elimination)
953 42 : _constraint_elimination->clear();
954 35602 : }
955 :
956 :
957 :
958 38844 : std::size_t DofMap::distribute_dofs (MeshBase & mesh)
959 : {
960 : // This function must be run on all processors at once
961 11218 : parallel_object_only();
962 :
963 : // Log how long it takes to distribute the degrees of freedom
964 33654 : LOG_SCOPE("distribute_dofs()", "DofMap");
965 :
966 11218 : libmesh_assert (mesh.is_prepared());
967 :
968 22436 : const processor_id_type proc_id = this->processor_id();
969 : #ifndef NDEBUG
970 11218 : const processor_id_type n_proc = this->n_processors();
971 : #endif
972 :
973 : // libmesh_assert_greater (this->n_variables(), 0);
974 11218 : libmesh_assert_less (proc_id, n_proc);
975 :
976 : // Data structure to ensure we can correctly combine
977 : // subdomain-restricted variables with constraining nodes from
978 : // different subdomains
979 : const std::map<const Node *, std::set<subdomain_id_type>>
980 : constraining_subdomains =
981 38846 : this->calculate_constraining_subdomains();
982 :
983 : // re-init in case the mesh has changed
984 38844 : this->reinit(mesh,
985 : constraining_subdomains);
986 :
987 : // By default distribute variables in a
988 : // var-major fashion, but allow run-time
989 : // specification
990 38837 : bool node_major_dofs = libMesh::on_command_line ("--node-major-dofs");
991 :
992 : // The DOF counter, will be incremented as we encounter
993 : // new degrees of freedom
994 38837 : dof_id_type next_free_dof = 0;
995 :
996 : // Clear the send list before we rebuild it
997 11216 : this->clear_send_list();
998 :
999 : // Set temporary DOF indices on this processor
1000 38837 : if (node_major_dofs)
1001 : this->distribute_local_dofs_node_major
1002 72 : (next_free_dof, mesh, constraining_subdomains);
1003 : else
1004 : this->distribute_local_dofs_var_major
1005 38765 : (next_free_dof, mesh, constraining_subdomains);
1006 :
1007 : // Get DOF counts on all processors
1008 38837 : const auto n_dofs = this->compute_dof_info(next_free_dof);
1009 :
1010 : // Clear all the current DOF indices
1011 : // (distribute_dofs expects them cleared!)
1012 38837 : this->invalidate_dofs(mesh);
1013 :
1014 38837 : next_free_dof = _first_df[proc_id];
1015 :
1016 : // Set permanent DOF indices on this processor
1017 38837 : if (node_major_dofs)
1018 : this->distribute_local_dofs_node_major
1019 72 : (next_free_dof, mesh, constraining_subdomains);
1020 : else
1021 : this->distribute_local_dofs_var_major
1022 38765 : (next_free_dof, mesh, constraining_subdomains);
1023 :
1024 11216 : libmesh_assert_equal_to (next_free_dof, _end_df[proc_id]);
1025 :
1026 : //------------------------------------------------------------
1027 : // At this point, all n_comp and dof_number values on local
1028 : // DofObjects should be correct, but a DistributedMesh might have
1029 : // incorrect values on non-local DofObjects. Let's request the
1030 : // correct values from each other processor.
1031 :
1032 50053 : if (this->n_processors() > 1)
1033 : {
1034 55856 : this->set_nonlocal_dof_objects(mesh.nodes_begin(),
1035 44752 : mesh.nodes_end(),
1036 : mesh, &DofMap::node_ptr);
1037 :
1038 55856 : this->set_nonlocal_dof_objects(mesh.elements_begin(),
1039 67077 : mesh.elements_end(),
1040 : mesh, &DofMap::elem_ptr);
1041 : }
1042 :
1043 : #ifdef DEBUG
1044 : {
1045 : const unsigned int
1046 11216 : sys_num = this->sys_number();
1047 :
1048 : // Processors should all agree on DoF ids for the newly numbered
1049 : // system.
1050 11216 : MeshTools::libmesh_assert_valid_dof_ids(mesh, sys_num);
1051 :
1052 : // DoF processor ids should match DofObject processor ids
1053 1950387 : for (auto & node : mesh.node_ptr_range())
1054 : {
1055 1939171 : DofObject const * const dofobj = node;
1056 1939171 : const processor_id_type obj_proc_id = dofobj->processor_id();
1057 :
1058 6669952 : for (auto v : make_range(dofobj->n_vars(sys_num)))
1059 7582950 : for (auto c : make_range(dofobj->n_comp(sys_num,v)))
1060 : {
1061 2852169 : const dof_id_type dofid = dofobj->dof_number(sys_num,v,c);
1062 2852169 : libmesh_assert_greater_equal (dofid, this->first_dof(obj_proc_id));
1063 2852169 : libmesh_assert_less (dofid, this->end_dof(obj_proc_id));
1064 : }
1065 : }
1066 :
1067 1058834 : for (auto & elem : mesh.element_ptr_range())
1068 : {
1069 1047618 : DofObject const * const dofobj = elem;
1070 1047618 : const processor_id_type obj_proc_id = dofobj->processor_id();
1071 :
1072 3261326 : for (auto v : make_range(dofobj->n_vars(sys_num)))
1073 3570906 : for (auto c : make_range(dofobj->n_comp(sys_num,v)))
1074 : {
1075 1357198 : const dof_id_type dofid = dofobj->dof_number(sys_num,v,c);
1076 1357198 : libmesh_assert_greater_equal (dofid, this->first_dof(obj_proc_id));
1077 1357198 : libmesh_assert_less (dofid, this->end_dof(obj_proc_id));
1078 : }
1079 : }
1080 : }
1081 : #endif
1082 :
1083 : // start finding SCALAR degrees of freedom
1084 : #ifdef LIBMESH_ENABLE_AMR
1085 38837 : _first_old_scalar_df = _first_scalar_df;
1086 : #endif
1087 11216 : _first_scalar_df.clear();
1088 38837 : _first_scalar_df.resize(this->n_variables(), DofObject::invalid_id);
1089 38837 : dof_id_type current_SCALAR_dof_index = n_dofs - n_SCALAR_dofs();
1090 :
1091 : // Calculate and cache the initial DoF indices for SCALAR variables.
1092 : // This is an O(N_vars) calculation so we want to do it once per
1093 : // renumbering rather than once per SCALAR_dof_indices() call
1094 :
1095 808715 : for (auto v : make_range(this->n_variables()))
1096 742257 : if (this->variable(v).type().family == SCALAR)
1097 : {
1098 141 : _first_scalar_df[v] = current_SCALAR_dof_index;
1099 141 : current_SCALAR_dof_index += this->variable(v).type().order.get_order();
1100 : }
1101 :
1102 : // Allow our GhostingFunctor objects to reinit if necessary
1103 77863 : for (const auto & gf : _algebraic_ghosting_functors)
1104 : {
1105 11270 : libmesh_assert(gf);
1106 39026 : gf->dofmap_reinit();
1107 : }
1108 :
1109 77674 : for (const auto & gf : _coupling_functors)
1110 : {
1111 11216 : libmesh_assert(gf);
1112 38837 : gf->dofmap_reinit();
1113 : }
1114 :
1115 : // Note that in the add_neighbors_to_send_list nodes on processor
1116 : // boundaries that are shared by multiple elements are added for
1117 : // each element.
1118 38837 : this->add_neighbors_to_send_list(mesh);
1119 :
1120 : // Here we used to clean up that data structure; now System and
1121 : // EquationSystems call that for us, after we've added constraint
1122 : // dependencies to the send_list too.
1123 : // this->sort_send_list ();
1124 :
1125 50053 : return n_dofs;
1126 : }
1127 :
1128 :
1129 : template <typename T, std::enable_if_t<std::is_same_v<T, dof_id_type> ||
1130 : std::is_same_v<T, std::vector<dof_id_type>>, int>>
1131 172 : void DofMap::local_variable_indices(T & idx,
1132 : const MeshBase & mesh,
1133 : unsigned int var_num) const
1134 : {
1135 : // Only used if T == dof_id_type to keep track of the greatest dof we've seen
1136 44 : dof_id_type greatest = 0;
1137 :
1138 : if constexpr (std::is_same_v<T, dof_id_type>)
1139 50 : idx = 0;
1140 : else if constexpr (std::is_same_v<T, std::vector<dof_id_type>>)
1141 40 : idx.clear();
1142 :
1143 : // Count dofs in the *exact* order that distribute_dofs numbered
1144 : // them, so that we can assume ascending indices and use push_back
1145 : // instead of find+insert.
1146 :
1147 88 : const unsigned int sys_num = this->sys_number();
1148 :
1149 : // If this isn't a SCALAR variable, we need to find all its field
1150 : // dofs on the mesh
1151 172 : if (this->variable_type(var_num).family != SCALAR)
1152 : {
1153 172 : const Variable & var(this->variable(var_num));
1154 :
1155 54586 : for (auto & elem : mesh.active_local_element_ptr_range())
1156 : {
1157 31677 : if (!var.active_on_subdomain(elem->subdomain_id()))
1158 48 : continue;
1159 :
1160 : // Only count dofs connected to active
1161 : // elements on this processor.
1162 31613 : const unsigned int n_nodes = elem->n_nodes();
1163 :
1164 : // First get any new nodal DOFS
1165 408677 : for (unsigned int n=0; n<n_nodes; n++)
1166 : {
1167 377064 : const Node & node = elem->node_ref(n);
1168 :
1169 417976 : if (node.processor_id() != this->processor_id())
1170 9349 : continue;
1171 :
1172 366919 : const unsigned int n_comp = node.n_comp(sys_num, var_num);
1173 544939 : for(unsigned int i=0; i<n_comp; i++)
1174 : {
1175 178020 : const dof_id_type index = node.dof_number(sys_num,var_num,i);
1176 32636 : libmesh_assert (this->local_index(index));
1177 :
1178 : if constexpr (std::is_same_v<T, dof_id_type>)
1179 : {
1180 80090 : if (idx == 0 || index > greatest)
1181 37242 : { idx++; greatest = index; }
1182 : }
1183 : else if constexpr (std::is_same_v<T, std::vector<dof_id_type>>)
1184 : {
1185 97930 : if (idx.empty() || index > idx.back())
1186 43509 : idx.push_back(index);
1187 : }
1188 : }
1189 : }
1190 :
1191 : // Next get any new element DOFS
1192 31613 : const unsigned int n_comp = elem->n_comp(sys_num, var_num);
1193 31613 : for (unsigned int i=0; i<n_comp; i++)
1194 : {
1195 0 : const dof_id_type index = elem->dof_number(sys_num,var_num,i);
1196 :
1197 : if constexpr (std::is_same_v<T, dof_id_type>)
1198 : {
1199 0 : if (idx == 0 || index > greatest)
1200 0 : { idx++; greatest = index; }
1201 : }
1202 : else if constexpr (std::is_same_v<T, std::vector<dof_id_type>>)
1203 : {
1204 0 : if (idx.empty() || index > idx.back())
1205 0 : idx.push_back(index);
1206 : }
1207 : }
1208 : } // done looping over elements
1209 :
1210 :
1211 : // we may have missed assigning DOFs to nodes that we own
1212 : // but to which we have no connected elements matching our
1213 : // variable restriction criterion. this will happen, for example,
1214 : // if variable V is restricted to subdomain S. We may not own
1215 : // any elements which live in S, but we may own nodes which are
1216 : // *connected* to elements which do. in this scenario these nodes
1217 : // will presently have unnumbered DOFs. we need to take care of
1218 : // them here since we own them and no other processor will touch them.
1219 334282 : for (const auto & node : mesh.local_node_ptr_range())
1220 : {
1221 43286 : libmesh_assert(node);
1222 :
1223 210235 : const unsigned int n_comp = node->n_comp(sys_num, var_num);
1224 290995 : for (unsigned int i=0; i<n_comp; i++)
1225 : {
1226 80760 : const dof_id_type index = node->dof_number(sys_num,var_num,i);
1227 :
1228 : if constexpr (std::is_same_v<T, dof_id_type>)
1229 : {
1230 37242 : if (idx == 0 || index > greatest)
1231 0 : { idx++; greatest = index; }
1232 : }
1233 : else if constexpr (std::is_same_v<T, std::vector<dof_id_type>>)
1234 : {
1235 43518 : if (idx.empty() || index > idx.back())
1236 9 : idx.push_back(index);
1237 : }
1238 : }
1239 : }
1240 : }
1241 : // Otherwise, count up the SCALAR dofs, if we're on the processor
1242 : // that holds this SCALAR variable
1243 0 : else if (this->processor_id() == (this->n_processors()-1))
1244 : {
1245 0 : std::vector<dof_id_type> di_scalar;
1246 0 : this->SCALAR_dof_indices(di_scalar,var_num);
1247 :
1248 : if constexpr (std::is_same_v<T, dof_id_type>)
1249 0 : idx += std::distance(di_scalar.begin(), di_scalar.end());
1250 : else if constexpr (std::is_same_v<T, std::vector<dof_id_type>>)
1251 0 : idx.insert(idx.end(), di_scalar.begin(), di_scalar.end());
1252 : }
1253 172 : }
1254 :
1255 : template void DofMap::local_variable_indices(dof_id_type &,
1256 : const MeshBase &,
1257 : unsigned int) const;
1258 :
1259 : template void DofMap::local_variable_indices(std::vector<dof_id_type> &,
1260 : const MeshBase &,
1261 : unsigned int) const;
1262 :
1263 :
1264 : std::map<const Node *, std::set<subdomain_id_type>>
1265 38844 : DofMap::calculate_constraining_subdomains()
1266 : {
1267 11218 : std::map<const Node *, std::set<subdomain_id_type>> constraining_subdomains;
1268 38844 : const auto & constraint_rows = _mesh.get_constraint_rows();
1269 :
1270 : // We can't just loop over constraint rows here because we need
1271 : // element subdomain ids for the constrained nodes, but we don't
1272 : // want an extra loop if there are no constraint rows.
1273 38844 : if (!constraint_rows.empty())
1274 51076 : for (auto & elem : _mesh.active_element_ptr_range())
1275 : {
1276 35697 : const subdomain_id_type sbdid = elem->subdomain_id();
1277 :
1278 189772 : for (const Node & node : elem->node_ref_range())
1279 : {
1280 154075 : if (auto it = constraint_rows.find(&node);
1281 44334 : it != constraint_rows.end())
1282 : {
1283 641832 : for (const auto & [pr, val] : it->second)
1284 : {
1285 : const Node * spline_node =
1286 507093 : pr.first->node_ptr(pr.second);
1287 :
1288 507093 : constraining_subdomains[spline_node].insert(sbdid);
1289 : }
1290 : }
1291 : }
1292 64 : }
1293 :
1294 38844 : return constraining_subdomains;
1295 : }
1296 :
1297 :
1298 144 : void DofMap::distribute_local_dofs_node_major
1299 : (dof_id_type & next_free_dof,
1300 : MeshBase & mesh,
1301 : const std::map<const Node *, std::set<subdomain_id_type>> &
1302 : constraining_subdomains)
1303 : {
1304 96 : const unsigned int sys_num = this->sys_number();
1305 48 : const unsigned int n_var_groups = this->n_variable_groups();
1306 :
1307 : // This is the common case and we want to optimize for it
1308 : const bool constraining_subdomains_empty =
1309 48 : constraining_subdomains.empty();
1310 :
1311 : // Our numbering here must be kept consistent with the numbering
1312 : // scheme assumed by DofMap::local_variable_indices!
1313 :
1314 : //-------------------------------------------------------------------------
1315 : // First count and assign temporary numbers to local dofs
1316 25296 : for (auto & elem : mesh.active_local_element_ptr_range())
1317 : {
1318 : // Only number dofs connected to active
1319 : // elements on this processor.
1320 18792 : const unsigned int n_nodes = elem->n_nodes();
1321 :
1322 18792 : const subdomain_id_type sbdid = elem->subdomain_id();
1323 :
1324 : // First number the nodal DOFS
1325 187920 : for (unsigned int n=0; n<n_nodes; n++)
1326 : {
1327 112752 : Node & node = elem->node_ref(n);
1328 :
1329 507384 : for (unsigned vg=0; vg<n_var_groups; vg++)
1330 : {
1331 112752 : const VariableGroup & vg_description(this->variable_group(vg));
1332 :
1333 338256 : if (vg_description.type().family == SCALAR)
1334 0 : continue;
1335 :
1336 : bool active_on_node =
1337 225504 : vg_description.active_on_subdomain(sbdid);
1338 :
1339 : // Are we at least active indirectly here?
1340 338256 : if (!active_on_node && !constraining_subdomains_empty)
1341 0 : if (auto it = constraining_subdomains.find(&node);
1342 0 : it != constraining_subdomains.end())
1343 0 : for (auto s : it->second)
1344 0 : if (vg_description.active_on_subdomain(s))
1345 : {
1346 0 : active_on_node = true;
1347 0 : break;
1348 : }
1349 :
1350 338256 : if (active_on_node)
1351 : {
1352 : // assign dof numbers (all at once) if this is
1353 : // our node and if they aren't already there
1354 275616 : if ((node.n_comp_group(sys_num,vg) > 0) &&
1355 417874 : (node.processor_id() == this->processor_id()) &&
1356 79618 : (node.vg_dof_base(sys_num,vg) ==
1357 : DofObject::invalid_id))
1358 : {
1359 100368 : node.set_vg_dof_base(sys_num, vg,
1360 : next_free_dof);
1361 100368 : next_free_dof += (vg_description.n_variables()*
1362 33456 : node.n_comp_group(sys_num,vg));
1363 : //node.debug_buffer();
1364 : }
1365 : }
1366 : }
1367 : }
1368 :
1369 : // Now number the element DOFS
1370 56376 : for (unsigned vg=0; vg<n_var_groups; vg++)
1371 : {
1372 12528 : const VariableGroup & vg_description(this->variable_group(vg));
1373 :
1374 50112 : if ((vg_description.type().family != SCALAR) &&
1375 37584 : (vg_description.active_on_subdomain(elem->subdomain_id())))
1376 37584 : if (elem->n_comp_group(sys_num,vg) > 0)
1377 : {
1378 0 : libmesh_assert_equal_to (elem->vg_dof_base(sys_num,vg),
1379 : DofObject::invalid_id);
1380 :
1381 0 : elem->set_vg_dof_base(sys_num,
1382 : vg,
1383 : next_free_dof);
1384 :
1385 0 : next_free_dof += (vg_description.n_variables()*
1386 0 : elem->n_comp_group(sys_num,vg));
1387 : }
1388 : }
1389 48 : } // done looping over elements
1390 :
1391 :
1392 : // we may have missed assigning DOFs to nodes that we own
1393 : // but to which we have no connected elements matching our
1394 : // variable restriction criterion. this will happen, for example,
1395 : // if variable V is restricted to subdomain S. We may not own
1396 : // any elements which live in S, but we may own nodes which are
1397 : // *connected* to elements which do. in this scenario these nodes
1398 : // will presently have unnumbered DOFs. we need to take care of
1399 : // them here since we own them and no other processor will touch them.
1400 198672 : for (auto & node : mesh.local_node_ptr_range())
1401 446472 : for (unsigned vg=0; vg<n_var_groups; vg++)
1402 : {
1403 99216 : const VariableGroup & vg_description(this->variable_group(vg));
1404 :
1405 396864 : if (node->n_comp_group(sys_num,vg))
1406 100368 : if (node->vg_dof_base(sys_num,vg) == DofObject::invalid_id)
1407 : {
1408 0 : node->set_vg_dof_base (sys_num,
1409 : vg,
1410 : next_free_dof);
1411 :
1412 0 : next_free_dof += (vg_description.n_variables()*
1413 0 : node->n_comp(sys_num,vg));
1414 : }
1415 48 : }
1416 :
1417 144 : this->distribute_scalar_dofs(next_free_dof);
1418 :
1419 : #ifdef DEBUG
1420 48 : this->assert_no_nodes_missed(mesh);
1421 : #endif // DEBUG
1422 144 : }
1423 :
1424 :
1425 :
1426 77530 : void DofMap::distribute_local_dofs_var_major
1427 : (dof_id_type & next_free_dof,
1428 : MeshBase & mesh,
1429 : const std::map<const Node *, std::set<subdomain_id_type>> &
1430 : constraining_subdomains)
1431 : {
1432 44768 : const unsigned int sys_num = this->sys_number();
1433 22384 : const unsigned int n_var_groups = this->n_variable_groups();
1434 :
1435 : // This is the common case and we want to optimize for it
1436 : const bool constraining_subdomains_empty =
1437 22384 : constraining_subdomains.empty();
1438 :
1439 : // Our numbering here must be kept consistent with the numbering
1440 : // scheme assumed by DofMap::local_variable_indices!
1441 :
1442 : //-------------------------------------------------------------------------
1443 : // First count and assign temporary numbers to local dofs
1444 157078 : for (unsigned vg=0; vg<n_var_groups; vg++)
1445 : {
1446 22952 : const VariableGroup & vg_description(this->variable_group(vg));
1447 :
1448 22952 : const unsigned int n_vars_in_group = vg_description.n_variables();
1449 :
1450 : // Skip the SCALAR dofs
1451 79548 : if (vg_description.type().family == SCALAR)
1452 198 : continue;
1453 :
1454 5032256 : for (auto & elem : mesh.active_local_element_ptr_range())
1455 : {
1456 : // Only number dofs connected to active elements on this
1457 : // processor and only for variables which are active on on
1458 : // this element's subdomain or which are active on the
1459 : // subdomain of a node constrained by this node.
1460 : const bool active_on_elem =
1461 3319914 : vg_description.active_on_subdomain(elem->subdomain_id());
1462 :
1463 : // If there's no way we're active on this element then we're
1464 : // done
1465 3319914 : if (!active_on_elem && constraining_subdomains_empty)
1466 5976 : continue;
1467 :
1468 3311898 : const unsigned int n_nodes = elem->n_nodes();
1469 :
1470 : // First number the nodal DOFS
1471 29383228 : for (unsigned int n=0; n<n_nodes; n++)
1472 : {
1473 26071330 : Node & node = elem->node_ref(n);
1474 :
1475 6509344 : bool active_on_node = active_on_elem;
1476 26071330 : if (!active_on_node)
1477 0 : if (auto it = constraining_subdomains.find(&node);
1478 0 : it != constraining_subdomains.end())
1479 0 : for (auto s : it->second)
1480 0 : if (vg_description.active_on_subdomain(s))
1481 : {
1482 0 : active_on_node = true;
1483 0 : break;
1484 : }
1485 :
1486 13018688 : if (!active_on_node)
1487 0 : continue;
1488 :
1489 : // assign dof numbers (all at once) if this is
1490 : // our node and if they aren't already there
1491 16409588 : if ((node.n_comp_group(sys_num,vg) > 0) &&
1492 29484008 : (node.processor_id() == this->processor_id()) &&
1493 3412678 : (node.vg_dof_base(sys_num,vg) ==
1494 : DofObject::invalid_id))
1495 : {
1496 4359314 : node.set_vg_dof_base(sys_num, vg, next_free_dof);
1497 :
1498 4359314 : next_free_dof += (n_vars_in_group*
1499 1138880 : node.n_comp_group(sys_num,vg));
1500 : }
1501 : }
1502 :
1503 : // Now number the element DOFS
1504 4181476 : if (elem->n_comp_group(sys_num,vg) > 0)
1505 : {
1506 235450 : libmesh_assert_equal_to (elem->vg_dof_base(sys_num,vg),
1507 : DofObject::invalid_id);
1508 :
1509 883954 : elem->set_vg_dof_base(sys_num,
1510 : vg,
1511 : next_free_dof);
1512 :
1513 883954 : next_free_dof += (n_vars_in_group*
1514 235450 : elem->n_comp_group(sys_num,vg));
1515 : }
1516 33530 : } // end loop on elements
1517 :
1518 : // we may have missed assigning DOFs to nodes that we own
1519 : // but to which we have no connected elements matching our
1520 : // variable restriction criterion. this will happen, for example,
1521 : // if variable V is restricted to subdomain S. We may not own
1522 : // any elements which live in S, but we may own nodes which are
1523 : // *connected* to elements which do. in this scenario these nodes
1524 : // will presently have unnumbered DOFs. we need to take care of
1525 : // them here since we own them and no other processor will touch them.
1526 13274128 : for (auto & node : mesh.local_node_ptr_range())
1527 10942032 : if (node->n_comp_group(sys_num,vg))
1528 4361360 : if (node->vg_dof_base(sys_num,vg) == DofObject::invalid_id)
1529 : {
1530 2046 : node->set_vg_dof_base (sys_num,
1531 : vg,
1532 : next_free_dof);
1533 :
1534 2046 : next_free_dof += (n_vars_in_group*
1535 682 : node->n_comp_group(sys_num,vg));
1536 33530 : }
1537 : } // end loop on variable groups
1538 :
1539 77530 : this->distribute_scalar_dofs(next_free_dof);
1540 :
1541 : #ifdef DEBUG
1542 22384 : this->assert_no_nodes_missed(mesh);
1543 : #endif
1544 77530 : }
1545 :
1546 :
1547 :
1548 77674 : void DofMap::distribute_scalar_dofs(dof_id_type & next_free_dof)
1549 : {
1550 77674 : this->_n_SCALAR_dofs = 0;
1551 157510 : for (auto vg : make_range(this->n_variable_groups()))
1552 : {
1553 23048 : const VariableGroup & vg_description(this->variable_group(vg));
1554 :
1555 79836 : if (vg_description.type().family == SCALAR)
1556 : {
1557 282 : this->_n_SCALAR_dofs += (vg_description.n_variables()*
1558 282 : vg_description.type().order.get_order());
1559 282 : continue;
1560 : }
1561 : }
1562 :
1563 : // Only increment next_free_dof if we're on the processor
1564 : // that holds this SCALAR variable
1565 100106 : if (this->processor_id() == (this->n_processors()-1))
1566 44138 : next_free_dof += _n_SCALAR_dofs;
1567 77674 : }
1568 :
1569 :
1570 :
1571 : #ifdef DEBUG
1572 22432 : void DofMap::assert_no_nodes_missed(MeshBase & mesh)
1573 : {
1574 22432 : MeshTools::libmesh_assert_valid_procids<Node>(mesh);
1575 :
1576 1972218 : for (auto & node : mesh.local_node_ptr_range())
1577 : {
1578 1949786 : unsigned int n_var_g = node->n_var_groups(this->sys_number());
1579 4359068 : for (unsigned int vg=0; vg != n_var_g; ++vg)
1580 : {
1581 : unsigned int n_comp_g =
1582 2409282 : node->n_comp_group(this->sys_number(), vg);
1583 2409282 : dof_id_type my_first_dof = n_comp_g ?
1584 1173018 : node->vg_dof_base(this->sys_number(), vg) : 0;
1585 2409282 : libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
1586 : }
1587 : }
1588 22432 : }
1589 : #endif // DEBUG
1590 :
1591 :
1592 : void
1593 1257455 : DofMap::
1594 : merge_ghost_functor_outputs(GhostingFunctor::map_type & elements_to_ghost,
1595 : CouplingMatricesSet & temporary_coupling_matrices,
1596 : const GhostingFunctorIterator & gf_begin,
1597 : const GhostingFunctorIterator & gf_end,
1598 : const MeshBase::const_element_iterator & elems_begin,
1599 : const MeshBase::const_element_iterator & elems_end,
1600 : processor_id_type p)
1601 : {
1602 2548956 : for (const auto & gf : as_range(gf_begin, gf_end))
1603 : {
1604 706948 : GhostingFunctor::map_type more_elements_to_ghost;
1605 :
1606 353474 : libmesh_assert(gf);
1607 1291501 : (*gf)(elems_begin, elems_end, p, more_elements_to_ghost);
1608 :
1609 : // A GhostingFunctor should only return active elements, but
1610 : // I forgot to *document* that, so let's go as easy as we
1611 : // can on functors that return inactive elements.
1612 : #if defined(LIBMESH_ENABLE_DEPRECATED) && defined(LIBMESH_ENABLE_AMR)
1613 706948 : std::vector<std::pair<const Elem*, const CouplingMatrix*>> children_to_couple;
1614 1653582 : for (auto it = more_elements_to_ghost.begin();
1615 3838645 : it != more_elements_to_ghost.end();)
1616 : {
1617 2547144 : const Elem * elem = it->first;
1618 715555 : if (!elem->active())
1619 : {
1620 : libmesh_deprecated();
1621 0 : std::vector<const Elem*> children_to_ghost;
1622 0 : elem->active_family_tree(children_to_ghost,
1623 : /*reset=*/ false);
1624 0 : for (const Elem * child : children_to_ghost)
1625 0 : if (child->processor_id() != p)
1626 0 : children_to_couple.emplace_back(child, it->second);
1627 :
1628 0 : it = more_elements_to_ghost.erase(it);
1629 : }
1630 : else
1631 715555 : ++it;
1632 : }
1633 353474 : more_elements_to_ghost.insert(children_to_couple.begin(),
1634 : children_to_couple.end());
1635 : #endif
1636 :
1637 3838645 : for (const auto & [elem, elem_cm] : more_elements_to_ghost)
1638 : {
1639 : // At this point we should only have active elements, even
1640 : // if we had to fix up gf output to get here.
1641 715555 : libmesh_assert(elem->active());
1642 :
1643 2547144 : if (const auto existing_it = elements_to_ghost.find(elem);
1644 715555 : existing_it == elements_to_ghost.end())
1645 1774961 : elements_to_ghost.emplace(elem, elem_cm);
1646 : else
1647 : {
1648 82581 : if (existing_it->second)
1649 : {
1650 0 : if (elem_cm)
1651 : {
1652 : // If this isn't already a temporary
1653 : // then we need to make one so we'll
1654 : // have a non-const matrix to merge
1655 0 : if (temporary_coupling_matrices.empty() ||
1656 0 : !temporary_coupling_matrices.count(existing_it->second))
1657 : {
1658 : // Make copy. This just calls the
1659 : // compiler-generated copy constructor
1660 : // because the CouplingMatrix class does not
1661 : // define a custom copy constructor.
1662 0 : auto result_pr = temporary_coupling_matrices.insert(std::make_unique<CouplingMatrix>(*existing_it->second));
1663 0 : existing_it->second = result_pr.first->get();
1664 : }
1665 :
1666 : // Merge elem_cm into existing CouplingMatrix
1667 0 : const_cast<CouplingMatrix &>(*existing_it->second) &= *elem_cm;
1668 : }
1669 : else // elem_cm == nullptr
1670 : {
1671 : // Any existing_it matrix merged with a full
1672 : // matrix (symbolized as nullptr) gives another
1673 : // full matrix (symbolizable as nullptr).
1674 :
1675 : // So if existing_it->second is a temporary then
1676 : // we don't need it anymore; we might as well
1677 : // remove it to keep the set of temporaries
1678 : // small.
1679 0 : if (const auto temp_it = temporary_coupling_matrices.find(existing_it->second);
1680 0 : temp_it != temporary_coupling_matrices.end())
1681 0 : temporary_coupling_matrices.erase(temp_it);
1682 :
1683 0 : existing_it->second = nullptr;
1684 : }
1685 : }
1686 : // else we have a nullptr already, then we have a full
1687 : // coupling matrix, already, and merging with anything
1688 : // else won't change that, so we're done.
1689 : }
1690 : }
1691 : }
1692 1257455 : }
1693 :
1694 :
1695 :
1696 38873 : void DofMap::add_neighbors_to_send_list(MeshBase & mesh)
1697 : {
1698 22456 : LOG_SCOPE("add_neighbors_to_send_list()", "DofMap");
1699 :
1700 : // Return immediately if there's no ghost data
1701 50101 : if (this->n_processors() == 1)
1702 5301 : return;
1703 :
1704 33572 : const unsigned int n_var = this->n_variables();
1705 :
1706 : MeshBase::const_element_iterator local_elem_it
1707 44800 : = mesh.active_local_elements_begin();
1708 : const MeshBase::const_element_iterator local_elem_end
1709 67144 : = mesh.active_local_elements_end();
1710 :
1711 22456 : GhostingFunctor::map_type elements_to_send;
1712 22456 : DofMap::CouplingMatricesSet temporary_coupling_matrices;
1713 :
1714 : // We need to add dofs to the send list if they've been directly
1715 : // requested by an algebraic ghosting functor or they've been
1716 : // indirectly requested by a coupling functor.
1717 44800 : this->merge_ghost_functor_outputs(elements_to_send,
1718 : temporary_coupling_matrices,
1719 44688 : this->algebraic_ghosting_functors_begin(),
1720 33572 : this->algebraic_ghosting_functors_end(),
1721 : local_elem_it, local_elem_end, mesh.processor_id());
1722 :
1723 44800 : this->merge_ghost_functor_outputs(elements_to_send,
1724 : temporary_coupling_matrices,
1725 44688 : this->coupling_functors_begin(),
1726 33572 : this->coupling_functors_end(),
1727 : local_elem_it, local_elem_end, mesh.processor_id());
1728 :
1729 : // Making a list of non-zero coupling matrix columns is an
1730 : // O(N_var^2) operation. We cache it so we only have to do it once
1731 : // per CouplingMatrix and not once per element.
1732 : std::map<const CouplingMatrix *, std::vector<unsigned int>>
1733 22456 : column_variable_lists;
1734 :
1735 164017 : for (const auto & [partner, ghost_coupling] : elements_to_send)
1736 : {
1737 : // We asked ghosting functors not to give us local elements
1738 44842 : libmesh_assert_not_equal_to
1739 : (partner->processor_id(), this->processor_id());
1740 :
1741 : // Loop over any present coupling matrix column variables if we
1742 : // have a coupling matrix, or just add all variables to
1743 : // send_list if not.
1744 130445 : if (ghost_coupling)
1745 : {
1746 6 : libmesh_assert_equal_to (ghost_coupling->size(), n_var);
1747 :
1748 : // Try to find a cached list of column variables.
1749 : std::map<const CouplingMatrix *, std::vector<unsigned int>>::const_iterator
1750 6 : column_variable_list = column_variable_lists.find(ghost_coupling);
1751 :
1752 : // If we didn't find it, then we need to create it.
1753 18 : if (column_variable_list == column_variable_lists.end())
1754 : {
1755 : auto inserted_variable_list_pair =
1756 12 : column_variable_lists.emplace(ghost_coupling, std::vector<unsigned int>());
1757 4 : column_variable_list = inserted_variable_list_pair.first;
1758 :
1759 : std::vector<unsigned int> & new_variable_list =
1760 12 : inserted_variable_list_pair.first->second;
1761 :
1762 16 : std::vector<unsigned char> has_variable(n_var, false);
1763 :
1764 60 : for (unsigned int vi = 0; vi != n_var; ++vi)
1765 : {
1766 48 : ConstCouplingRow ccr(vi, *ghost_coupling);
1767 :
1768 108 : for (const auto & vj : ccr)
1769 80 : has_variable[vj] = true;
1770 : }
1771 60 : for (unsigned int vj = 0; vj != n_var; ++vj)
1772 : {
1773 64 : if (has_variable[vj])
1774 48 : new_variable_list.push_back(vj);
1775 : }
1776 : }
1777 :
1778 : const std::vector<unsigned int> & variable_list =
1779 6 : column_variable_list->second;
1780 :
1781 90 : for (const auto & vj : variable_list)
1782 : {
1783 48 : std::vector<dof_id_type> di;
1784 72 : this->dof_indices (partner, di, vj);
1785 :
1786 : // Insert the remote DOF indices into the send list
1787 216 : for (auto d : di)
1788 192 : if (d != DofObject::invalid_id &&
1789 96 : !this->local_index(d))
1790 : {
1791 48 : libmesh_assert_less(d, this->n_dofs());
1792 144 : _send_list.push_back(d);
1793 : }
1794 : }
1795 : }
1796 : else
1797 : {
1798 89672 : std::vector<dof_id_type> di;
1799 130427 : this->dof_indices (partner, di);
1800 :
1801 : // Insert the remote DOF indices into the send list
1802 2811559 : for (const auto & dof : di)
1803 3580169 : if (dof != DofObject::invalid_id &&
1804 1782095 : !this->local_index(dof))
1805 : {
1806 728515 : libmesh_assert_less(dof, this->n_dofs());
1807 2171788 : _send_list.push_back(dof);
1808 : }
1809 : }
1810 :
1811 : }
1812 :
1813 : // We're now done with any merged coupling matrices we had to create.
1814 11228 : temporary_coupling_matrices.clear();
1815 :
1816 : //-------------------------------------------------------------------------
1817 : // Our coupling functors added dofs from neighboring elements to the
1818 : // send list, but we may still need to add non-local dofs from local
1819 : // elements.
1820 : //-------------------------------------------------------------------------
1821 :
1822 : // Loop over the active local elements, adding all active elements
1823 : // that neighbor an active local element to the send list.
1824 1654596 : for ( ; local_elem_it != local_elem_end; ++local_elem_it)
1825 : {
1826 1217456 : const Elem * elem = *local_elem_it;
1827 :
1828 813888 : std::vector<dof_id_type> di;
1829 1217456 : this->dof_indices (elem, di);
1830 :
1831 : // Insert the remote DOF indices into the send list
1832 13714580 : for (const auto & dof : di)
1833 16660845 : if (dof != DofObject::invalid_id &&
1834 8333403 : !this->local_index(dof))
1835 : {
1836 188089 : libmesh_assert_less(dof, this->n_dofs());
1837 563722 : _send_list.push_back(dof);
1838 : }
1839 : }
1840 : }
1841 :
1842 :
1843 :
1844 38860 : void DofMap::prepare_send_list ()
1845 : {
1846 22448 : LOG_SCOPE("prepare_send_list()", "DofMap");
1847 :
1848 : // Return immediately if there's no ghost data
1849 50084 : if (this->n_processors() == 1)
1850 5300 : return;
1851 :
1852 : // Check to see if we have any extra stuff to add to the send_list
1853 33560 : if (_extra_send_list_function)
1854 : {
1855 0 : if (_augment_send_list)
1856 : {
1857 0 : libmesh_here();
1858 0 : libMesh::out << "WARNING: You have specified both an extra send list function and object.\n"
1859 0 : << " Are you sure this is what you meant to do??"
1860 0 : << std::endl;
1861 : }
1862 :
1863 0 : _extra_send_list_function(_send_list, _extra_send_list_context);
1864 : }
1865 :
1866 33560 : if (_augment_send_list)
1867 0 : _augment_send_list->augment_send_list (_send_list);
1868 :
1869 : // First sort the send list. After this
1870 : // duplicated elements will be adjacent in the
1871 : // vector
1872 33560 : std::sort(_send_list.begin(), _send_list.end());
1873 :
1874 : // Now use std::unique to remove duplicate entries
1875 : std::vector<dof_id_type>::iterator new_end =
1876 33560 : std::unique (_send_list.begin(), _send_list.end());
1877 :
1878 : // Remove the end of the send_list. Use the "swap trick"
1879 : // from Effective STL
1880 55896 : std::vector<dof_id_type> (_send_list.begin(), new_end).swap (_send_list);
1881 :
1882 : // Make sure the send list has nothing invalid in it.
1883 11224 : libmesh_assert(_send_list.empty() || _send_list.back() < this->n_dofs());
1884 : }
1885 :
1886 36 : void DofMap::reinit_send_list (MeshBase & mesh)
1887 : {
1888 12 : this->clear_send_list();
1889 36 : this->add_neighbors_to_send_list(mesh);
1890 :
1891 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
1892 : // This is assuming that we only need to recommunicate
1893 : // the constraints and no new ones have been added since
1894 : // a previous call to reinit_constraints.
1895 36 : this->process_constraints(mesh);
1896 : #endif
1897 36 : this->prepare_send_list();
1898 36 : }
1899 :
1900 14 : void DofMap::set_implicit_neighbor_dofs(bool implicit_neighbor_dofs)
1901 : {
1902 14 : _implicit_neighbor_dofs_initialized = true;
1903 14 : _implicit_neighbor_dofs = implicit_neighbor_dofs;
1904 14 : }
1905 :
1906 0 : void DofMap::set_verify_dirichlet_bc_consistency(bool val)
1907 : {
1908 0 : _verify_dirichlet_bc_consistency = val;
1909 0 : }
1910 :
1911 :
1912 79126 : bool DofMap::use_coupled_neighbor_dofs(const MeshBase & /*mesh*/) const
1913 : {
1914 : // If we were asked on the command line, then we need to
1915 : // include sensitivities between neighbor degrees of freedom
1916 : bool implicit_neighbor_dofs =
1917 79126 : libMesh::on_command_line ("--implicit-neighbor-dofs");
1918 :
1919 : // If the user specifies --implicit-neighbor-dofs 0, then
1920 : // presumably he knows what he is doing and we won't try to
1921 : // automatically turn it on even when all the variables are
1922 : // discontinuous.
1923 79126 : if (implicit_neighbor_dofs)
1924 : {
1925 : // No flag provided defaults to 'true'
1926 0 : int flag = 1;
1927 0 : flag = libMesh::command_line_next ("--implicit-neighbor-dofs", flag);
1928 :
1929 0 : if (!flag)
1930 : {
1931 : // The user said --implicit-neighbor-dofs 0, so he knows
1932 : // what he is doing and really doesn't want it.
1933 0 : return false;
1934 : }
1935 : }
1936 :
1937 : // Possibly override the commandline option, if set_implicit_neighbor_dofs
1938 : // has been called.
1939 79126 : if (_implicit_neighbor_dofs_initialized)
1940 : {
1941 42 : implicit_neighbor_dofs = _implicit_neighbor_dofs;
1942 :
1943 : // Again, if the user explicitly says implicit_neighbor_dofs = false,
1944 : // then we return here.
1945 42 : if (!implicit_neighbor_dofs)
1946 0 : return false;
1947 : }
1948 :
1949 : // Look at all the variables in this system. If every one is
1950 : // discontinuous then the user must be doing DG/FVM, so be nice
1951 : // and force implicit_neighbor_dofs=true.
1952 : {
1953 22852 : bool all_discontinuous_dofs = true;
1954 :
1955 : // We may call this method even without ever having initialized our data
1956 1566935 : for (auto var : index_range(this->_variables))
1957 1487809 : if (FEInterface::get_continuity(this->variable_type(var)) != DISCONTINUOUS)
1958 412894 : all_discontinuous_dofs = false;
1959 :
1960 79126 : if (all_discontinuous_dofs)
1961 11648 : implicit_neighbor_dofs = true;
1962 : }
1963 :
1964 22852 : return implicit_neighbor_dofs;
1965 : }
1966 :
1967 :
1968 :
1969 4206 : void DofMap::compute_sparsity(const MeshBase & mesh)
1970 : {
1971 7148 : _sp = this->build_sparsity(mesh, this->_constrained_sparsity_construction);
1972 :
1973 : // It is possible that some \p SparseMatrix implementations want to
1974 : // see the sparsity pattern before we throw it away. If so, we
1975 : // share a view of its arrays, and we pass it in to the matrices.
1976 8484 : for (const auto & mat : _matrices)
1977 : {
1978 5564 : mat->attach_sparsity_pattern (*_sp);
1979 4278 : if (_need_full_sparsity_pattern)
1980 0 : mat->update_sparsity_pattern (_sp->get_sparsity_pattern());
1981 : }
1982 : // If we don't need the full sparsity pattern anymore, free the
1983 : // parts of it we don't need.
1984 4206 : if (!_need_full_sparsity_pattern)
1985 4206 : _sp->clear_full_sparsity();
1986 4206 : }
1987 :
1988 :
1989 :
1990 37845 : void DofMap::clear_sparsity()
1991 : {
1992 10908 : _sp.reset();
1993 37845 : }
1994 :
1995 :
1996 :
1997 35 : void DofMap::remove_default_ghosting()
1998 : {
1999 35 : this->remove_coupling_functor(this->default_coupling());
2000 35 : this->remove_algebraic_ghosting_functor(this->default_algebraic_ghosting());
2001 35 : }
2002 :
2003 :
2004 :
2005 7 : void DofMap::add_default_ghosting()
2006 : {
2007 7 : this->add_coupling_functor(this->default_coupling());
2008 7 : this->add_algebraic_ghosting_functor(this->default_algebraic_ghosting());
2009 7 : }
2010 :
2011 :
2012 :
2013 : void
2014 71184 : DofMap::add_coupling_functor(GhostingFunctor & coupling_functor,
2015 : bool to_mesh)
2016 : {
2017 : // We used to implicitly support duplicate inserts to std::set
2018 : #ifdef LIBMESH_ENABLE_DEPRECATED
2019 : _coupling_functors.erase
2020 30320 : (std::remove(_coupling_functors.begin(),
2021 : _coupling_functors.end(),
2022 71184 : &coupling_functor),
2023 61296 : _coupling_functors.end());
2024 : #endif
2025 :
2026 : // We shouldn't have two copies of the same functor
2027 20432 : libmesh_assert(std::find(_coupling_functors.begin(),
2028 : _coupling_functors.end(),
2029 : &coupling_functor) ==
2030 : _coupling_functors.end());
2031 :
2032 71184 : _coupling_functors.push_back(&coupling_functor);
2033 71184 : coupling_functor.set_mesh(&_mesh);
2034 71184 : if (to_mesh)
2035 71177 : _mesh.add_ghosting_functor(coupling_functor);
2036 71184 : }
2037 :
2038 :
2039 :
2040 : void
2041 63 : DofMap::remove_coupling_functor(GhostingFunctor & coupling_functor)
2042 : {
2043 27 : auto raw_it = std::find(_coupling_functors.begin(),
2044 81 : _coupling_functors.end(), &coupling_functor);
2045 :
2046 : #ifndef LIBMESH_ENABLE_DEPRECATED
2047 : // We shouldn't be trying to remove a functor that isn't there
2048 : libmesh_assert(raw_it != _coupling_functors.end());
2049 : #else
2050 : // Our old API supported trying to remove a functor that isn't there
2051 63 : if (raw_it != _coupling_functors.end())
2052 : #endif
2053 63 : _coupling_functors.erase(raw_it);
2054 :
2055 : // We shouldn't have had two copies of the same functor
2056 18 : libmesh_assert(std::find(_coupling_functors.begin(),
2057 : _coupling_functors.end(),
2058 : &coupling_functor) ==
2059 : _coupling_functors.end());
2060 :
2061 63 : _mesh.remove_ghosting_functor(coupling_functor);
2062 :
2063 63 : if (const auto it = _shared_functors.find(&coupling_functor);
2064 18 : it != _shared_functors.end())
2065 0 : _shared_functors.erase(it);
2066 63 : }
2067 :
2068 :
2069 :
2070 : void
2071 71212 : DofMap::add_algebraic_ghosting_functor(GhostingFunctor & evaluable_functor,
2072 : bool to_mesh)
2073 : {
2074 : // We used to implicitly support duplicate inserts to std::set
2075 : #ifdef LIBMESH_ENABLE_DEPRECATED
2076 : _algebraic_ghosting_functors.erase
2077 30332 : (std::remove(_algebraic_ghosting_functors.begin(),
2078 : _algebraic_ghosting_functors.end(),
2079 71212 : &evaluable_functor),
2080 61320 : _algebraic_ghosting_functors.end());
2081 : #endif
2082 :
2083 : // We shouldn't have two copies of the same functor
2084 20440 : libmesh_assert(std::find(_algebraic_ghosting_functors.begin(),
2085 : _algebraic_ghosting_functors.end(),
2086 : &evaluable_functor) ==
2087 : _algebraic_ghosting_functors.end());
2088 :
2089 71212 : _algebraic_ghosting_functors.push_back(&evaluable_functor);
2090 71212 : evaluable_functor.set_mesh(&_mesh);
2091 71212 : if (to_mesh)
2092 71212 : _mesh.add_ghosting_functor(evaluable_functor);
2093 71212 : }
2094 :
2095 :
2096 :
2097 : void
2098 63 : DofMap::remove_algebraic_ghosting_functor(GhostingFunctor & evaluable_functor)
2099 : {
2100 27 : auto raw_it = std::find(_algebraic_ghosting_functors.begin(),
2101 : _algebraic_ghosting_functors.end(),
2102 81 : &evaluable_functor);
2103 :
2104 : #ifndef LIBMESH_ENABLE_DEPRECATED
2105 : // We shouldn't be trying to remove a functor that isn't there
2106 : libmesh_assert(raw_it != _algebraic_ghosting_functors.end());
2107 : #else
2108 : // Our old API supported trying to remove a functor that isn't there
2109 63 : if (raw_it != _algebraic_ghosting_functors.end())
2110 : #endif
2111 63 : _algebraic_ghosting_functors.erase(raw_it);
2112 :
2113 : // We shouldn't have had two copies of the same functor
2114 18 : libmesh_assert(std::find(_algebraic_ghosting_functors.begin(),
2115 : _algebraic_ghosting_functors.end(),
2116 : &evaluable_functor) ==
2117 : _algebraic_ghosting_functors.end());
2118 :
2119 63 : _mesh.remove_ghosting_functor(evaluable_functor);
2120 :
2121 63 : if (const auto it = _shared_functors.find(&evaluable_functor);
2122 18 : it != _shared_functors.end())
2123 0 : _shared_functors.erase(it);
2124 63 : }
2125 :
2126 :
2127 :
2128 0 : void DofMap::extract_local_vector (const NumericVector<Number> & Ug,
2129 : const std::vector<dof_id_type> & dof_indices_in,
2130 : DenseVectorBase<Number> & Ue) const
2131 : {
2132 0 : const unsigned int n_original_dofs = dof_indices_in.size();
2133 :
2134 : #ifdef LIBMESH_ENABLE_AMR
2135 :
2136 : // Trivial mapping
2137 0 : libmesh_assert_equal_to (dof_indices_in.size(), Ue.size());
2138 0 : bool has_constrained_dofs = false;
2139 :
2140 0 : for (unsigned int il=0; il != n_original_dofs; ++il)
2141 : {
2142 0 : const dof_id_type ig = dof_indices_in[il];
2143 :
2144 0 : if (this->is_constrained_dof (ig)) has_constrained_dofs = true;
2145 :
2146 0 : libmesh_assert_less (ig, Ug.size());
2147 :
2148 0 : Ue.el(il) = Ug(ig);
2149 : }
2150 :
2151 : // If the element has any constrained DOFs then we need
2152 : // to account for them in the mapping. This will handle
2153 : // the case that the input vector is not constrained.
2154 0 : if (has_constrained_dofs)
2155 : {
2156 : // Copy the input DOF indices.
2157 0 : std::vector<dof_id_type> constrained_dof_indices(dof_indices_in);
2158 :
2159 0 : DenseMatrix<Number> C;
2160 0 : DenseVector<Number> H;
2161 :
2162 0 : this->build_constraint_matrix_and_vector (C, H, constrained_dof_indices);
2163 :
2164 0 : libmesh_assert_equal_to (dof_indices_in.size(), C.m());
2165 0 : libmesh_assert_equal_to (constrained_dof_indices.size(), C.n());
2166 :
2167 : // zero-out Ue
2168 0 : Ue.zero();
2169 :
2170 : // compute Ue = C Ug, with proper mapping.
2171 0 : for (unsigned int i=0; i != n_original_dofs; i++)
2172 : {
2173 0 : Ue.el(i) = H(i);
2174 :
2175 : const unsigned int n_constrained =
2176 0 : cast_int<unsigned int>(constrained_dof_indices.size());
2177 0 : for (unsigned int j=0; j<n_constrained; j++)
2178 : {
2179 0 : const dof_id_type jg = constrained_dof_indices[j];
2180 :
2181 : // If Ug is a serial or ghosted vector, then this assert is
2182 : // overzealous. If Ug is a parallel vector, then this assert
2183 : // is redundant.
2184 : // libmesh_assert ((jg >= Ug.first_local_index()) &&
2185 : // (jg < Ug.last_local_index()));
2186 :
2187 0 : Ue.el(i) += C(i,j)*Ug(jg);
2188 : }
2189 : }
2190 0 : }
2191 :
2192 : #else
2193 :
2194 : // Trivial mapping
2195 :
2196 : libmesh_assert_equal_to (n_original_dofs, Ue.size());
2197 :
2198 : for (unsigned int il=0; il<n_original_dofs; il++)
2199 : {
2200 : const dof_id_type ig = dof_indices_in[il];
2201 :
2202 : libmesh_assert ((ig >= Ug.first_local_index()) && (ig < Ug.last_local_index()));
2203 :
2204 : Ue.el(il) = Ug(ig);
2205 : }
2206 :
2207 : #endif
2208 0 : }
2209 :
2210 37774333 : void DofMap::dof_indices (const Elem * const elem,
2211 : std::vector<dof_id_type> & di) const
2212 : {
2213 : // We now allow elem==nullptr to request just SCALAR dofs
2214 : // libmesh_assert(elem);
2215 :
2216 : // If we are asking for current indices on an element, it ought to
2217 : // be an active element (or a temporary side, which also thinks it's
2218 : // active)
2219 10923734 : libmesh_assert(!elem || elem->active());
2220 :
2221 : // dof_indices() is a relatively light-weight function that is
2222 : // called millions of times in normal codes. Therefore, it is not a
2223 : // good candidate for logging, since the cost of the logging code
2224 : // itself is roughly on par with the time required to call
2225 : // dof_indices().
2226 : // LOG_SCOPE("dof_indices()", "DofMap");
2227 :
2228 : // Clear the DOF indices vector
2229 10923734 : di.clear();
2230 :
2231 10923734 : const unsigned int n_var_groups = this->n_variable_groups();
2232 :
2233 : #ifdef DEBUG
2234 : // Check that sizes match in DEBUG mode
2235 10923734 : std::size_t tot_size = 0;
2236 : #endif
2237 :
2238 37774333 : if (elem && elem->type() == TRI3SUBDIVISION)
2239 : {
2240 : // Subdivision surface FE require the 1-ring around elem
2241 616 : const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
2242 :
2243 : // Ghost subdivision elements have no real dofs
2244 2104 : if (!sd_elem->is_ghost())
2245 : {
2246 : // Determine the nodes contributing to element elem
2247 1088 : std::vector<const Node *> elem_nodes;
2248 1888 : MeshTools::Subdivision::find_one_ring(sd_elem, elem_nodes);
2249 :
2250 : // Get the dof numbers
2251 3776 : for (unsigned int vg=0; vg<n_var_groups; vg++)
2252 : {
2253 544 : const VariableGroup & var = this->variable_group(vg);
2254 544 : const unsigned int vars_in_group = var.n_variables();
2255 :
2256 1888 : if (var.type().family == SCALAR &&
2257 0 : var.active_on_subdomain(elem->subdomain_id()))
2258 : {
2259 0 : for (unsigned int vig=0; vig != vars_in_group; ++vig)
2260 : {
2261 : #ifdef DEBUG
2262 0 : tot_size += var.type().order;
2263 : #endif
2264 0 : std::vector<dof_id_type> di_new;
2265 0 : this->SCALAR_dof_indices(di_new,var.number(vig));
2266 0 : di.insert( di.end(), di_new.begin(), di_new.end());
2267 : }
2268 : }
2269 : else
2270 7552 : for (unsigned int vig=0; vig != vars_in_group; ++vig)
2271 : {
2272 10560 : _dof_indices(*elem, elem->p_level(), di, vg, vig,
2273 1632 : elem_nodes.data(),
2274 : cast_int<unsigned int>(elem_nodes.size()),
2275 : var.number(vig)
2276 : #ifdef DEBUG
2277 : , tot_size
2278 : #endif
2279 : );
2280 : }
2281 : }
2282 : }
2283 :
2284 2104 : return;
2285 : }
2286 :
2287 : // Get the dof numbers for each variable
2288 37772229 : const unsigned int n_nodes = elem ? elem->n_nodes() : 0;
2289 77672374 : for (unsigned int vg=0; vg<n_var_groups; vg++)
2290 : {
2291 11479109 : const VariableGroup & var = this->variable_group(vg);
2292 11479109 : const unsigned int vars_in_group = var.n_variables();
2293 :
2294 39979531 : if (var.type().family == SCALAR &&
2295 212829 : (!elem ||
2296 292215 : var.active_on_subdomain(elem->subdomain_id())))
2297 : {
2298 584430 : for (unsigned int vig=0; vig != vars_in_group; ++vig)
2299 : {
2300 : #ifdef DEBUG
2301 79386 : tot_size += var.type().order;
2302 : #endif
2303 79386 : std::vector<dof_id_type> di_new;
2304 371601 : this->SCALAR_dof_indices(di_new,var.number(vig));
2305 292215 : di.insert( di.end(), di_new.begin(), di_new.end());
2306 : }
2307 : }
2308 39607930 : else if (elem)
2309 84861611 : for (unsigned int vig=0; vig != vars_in_group; ++vig)
2310 : {
2311 58152030 : _dof_indices(*elem, elem->p_level(), di, vg, vig,
2312 : elem->get_nodes(), n_nodes, var.number(vig)
2313 : #ifdef DEBUG
2314 : , tot_size
2315 : #endif
2316 : );
2317 : }
2318 : }
2319 :
2320 : #ifdef DEBUG
2321 10923118 : libmesh_assert_equal_to (tot_size, di.size());
2322 : #endif
2323 : }
2324 :
2325 :
2326 36876467 : void DofMap::dof_indices (const Elem * const elem,
2327 : std::vector<dof_id_type> & di,
2328 : const unsigned int vn,
2329 : int p_level) const
2330 : {
2331 36876467 : dof_indices(
2332 : elem,
2333 : di,
2334 : vn,
2335 182298 : [](const Elem &,
2336 : std::vector<dof_id_type> & dof_indices,
2337 195846 : const std::vector<dof_id_type> & scalar_dof_indices) {
2338 476067 : dof_indices.insert(dof_indices.end(), scalar_dof_indices.begin(), scalar_dof_indices.end());
2339 378144 : },
2340 : [](const Elem &,
2341 : unsigned int,
2342 : unsigned int,
2343 : std::vector<dof_id_type> & dof_indices,
2344 143537568 : const dof_id_type dof) { dof_indices.push_back(dof); },
2345 : p_level);
2346 36876467 : }
2347 :
2348 56 : void DofMap::array_dof_indices(const Elem * const elem,
2349 : std::vector<dof_id_type> & di,
2350 : const unsigned int vn,
2351 : int p_level) const
2352 : {
2353 30 : auto dof_indices_functor = [elem, p_level, this](std::vector<dof_id_type> & functor_di,
2354 110 : const unsigned int functor_vn) {
2355 70 : this->dof_indices(elem, functor_di, functor_vn, p_level);
2356 106 : };
2357 56 : this->array_dof_indices(dof_indices_functor, di, vn);
2358 56 : }
2359 :
2360 0 : void DofMap::array_dof_indices(const Node * const node,
2361 : std::vector<dof_id_type> & di,
2362 : const unsigned int vn) const
2363 : {
2364 : auto dof_indices_functor = [node, this](std::vector<dof_id_type> & functor_di,
2365 0 : const unsigned int functor_vn) {
2366 0 : this->dof_indices(node, functor_di, functor_vn);
2367 0 : };
2368 0 : this->array_dof_indices(dof_indices_functor, di, vn);
2369 0 : }
2370 :
2371 45684 : void DofMap::dof_indices (const Node * const node,
2372 : std::vector<dof_id_type> & di) const
2373 : {
2374 : // We allow node==nullptr to request just SCALAR dofs
2375 : // libmesh_assert(elem);
2376 :
2377 : // dof_indices() is a relatively light-weight function that is
2378 : // called millions of times in normal codes. Therefore, it is not a
2379 : // good candidate for logging, since the cost of the logging code
2380 : // itself is roughly on par with the time required to call
2381 : // dof_indices().
2382 : // LOG_SCOPE("dof_indices(Node)", "DofMap");
2383 :
2384 : // Clear the DOF indices vector
2385 15228 : di.clear();
2386 :
2387 15228 : const unsigned int n_var_groups = this->n_variable_groups();
2388 30456 : const unsigned int sys_num = this->sys_number();
2389 :
2390 : // Get the dof numbers
2391 137052 : for (unsigned int vg=0; vg<n_var_groups; vg++)
2392 : {
2393 30456 : const VariableGroup & var = this->variable_group(vg);
2394 30456 : const unsigned int vars_in_group = var.n_variables();
2395 :
2396 91368 : if (var.type().family == SCALAR)
2397 : {
2398 0 : for (unsigned int vig=0; vig != vars_in_group; ++vig)
2399 : {
2400 0 : std::vector<dof_id_type> di_new;
2401 0 : this->SCALAR_dof_indices(di_new,var.number(vig));
2402 0 : di.insert( di.end(), di_new.begin(), di_new.end());
2403 : }
2404 : }
2405 : else
2406 : {
2407 91368 : const int n_comp = node->n_comp_group(sys_num,vg);
2408 228420 : for (unsigned int vig=0; vig != vars_in_group; ++vig)
2409 : {
2410 248724 : for (int i=0; i != n_comp; ++i)
2411 : {
2412 : const dof_id_type d =
2413 111672 : node->dof_number(sys_num, vg, vig, i, n_comp);
2414 37224 : libmesh_assert_not_equal_to
2415 : (d, DofObject::invalid_id);
2416 111672 : di.push_back(d);
2417 : }
2418 : }
2419 : }
2420 : }
2421 45684 : }
2422 :
2423 :
2424 200 : void DofMap::dof_indices (const Node * const node,
2425 : std::vector<dof_id_type> & di,
2426 : const unsigned int vn) const
2427 : {
2428 200 : if (vn == libMesh::invalid_uint)
2429 : {
2430 0 : this->dof_indices(node, di);
2431 0 : return;
2432 : }
2433 :
2434 : // We allow node==nullptr to request just SCALAR dofs
2435 : // libmesh_assert(elem);
2436 :
2437 : // dof_indices() is a relatively light-weight function that is
2438 : // called millions of times in normal codes. Therefore, it is not a
2439 : // good candidate for logging, since the cost of the logging code
2440 : // itself is roughly on par with the time required to call
2441 : // dof_indices().
2442 : // LOG_SCOPE("dof_indices(Node)", "DofMap");
2443 :
2444 : // Clear the DOF indices vector
2445 50 : di.clear();
2446 :
2447 100 : const unsigned int sys_num = this->sys_number();
2448 :
2449 : // Get the dof numbers
2450 250 : const unsigned int vg = this->_variable_group_numbers[vn];
2451 50 : const VariableGroup & var = this->variable_group(vg);
2452 :
2453 200 : if (var.type().family == SCALAR)
2454 : {
2455 0 : std::vector<dof_id_type> di_new;
2456 0 : this->SCALAR_dof_indices(di_new,vn);
2457 0 : di.insert( di.end(), di_new.begin(), di_new.end());
2458 : }
2459 : else
2460 : {
2461 200 : const unsigned int vig = vn - var.number();
2462 200 : const int n_comp = node->n_comp_group(sys_num,vg);
2463 320 : for (int i=0; i != n_comp; ++i)
2464 : {
2465 : const dof_id_type d =
2466 120 : node->dof_number(sys_num, vg, vig, i, n_comp);
2467 30 : libmesh_assert_not_equal_to
2468 : (d, DofObject::invalid_id);
2469 120 : di.push_back(d);
2470 : }
2471 : }
2472 : }
2473 :
2474 :
2475 1627164 : void DofMap::dof_indices (const Elem & elem,
2476 : unsigned int n,
2477 : std::vector<dof_id_type> & di,
2478 : const unsigned int vn) const
2479 : {
2480 1627164 : this->_node_dof_indices(elem, n, elem.node_ref(n), di, vn);
2481 1627164 : }
2482 :
2483 :
2484 :
2485 : #ifdef LIBMESH_ENABLE_AMR
2486 :
2487 1515752 : void DofMap::old_dof_indices (const Elem & elem,
2488 : unsigned int n,
2489 : std::vector<dof_id_type> & di,
2490 : const unsigned int vn) const
2491 : {
2492 452311 : const DofObject & old_obj = elem.node_ref(n).get_old_dof_object_ref();
2493 1515752 : this->_node_dof_indices(elem, n, old_obj, di, vn);
2494 1515752 : }
2495 :
2496 : #endif // LIBMESH_ENABLE_AMR
2497 :
2498 :
2499 :
2500 3142916 : void DofMap::_node_dof_indices (const Elem & elem,
2501 : unsigned int n,
2502 : const DofObject & obj,
2503 : std::vector<dof_id_type> & di,
2504 : const unsigned int vn) const
2505 : {
2506 : // Half of this is a cut and paste of _dof_indices code below, but
2507 : // duplication actually seems cleaner than creating a helper
2508 : // function with a million arguments and hoping the compiler inlines
2509 : // it properly into one of our most highly trafficked functions.
2510 :
2511 : // dof_indices() is a relatively light-weight function; the cost of
2512 : // the logging code itself is roughly on par with the time required
2513 : // to call dof_indices().
2514 : // LOG_SCOPE("_node_dof_indices()", "DofMap");
2515 :
2516 1889230 : const unsigned int sys_num = this->sys_number();
2517 944615 : const auto [vg, vig] =
2518 2198301 : obj.var_to_vg_and_offset(sys_num,vn);
2519 2198301 : const unsigned int n_comp = obj.n_comp_group(sys_num,vg);
2520 :
2521 944615 : const VariableGroup & var = this->variable_group(vg);
2522 3142916 : FEType fe_type = var.type();
2523 : const bool extra_hanging_dofs =
2524 3142916 : FEInterface::extra_hanging_dofs(fe_type);
2525 :
2526 3142916 : const bool add_p_level = fe_type.p_refinement;
2527 :
2528 : // There is a potential problem with h refinement. Imagine a
2529 : // quad9 that has a linear FE on it. Then, on the hanging side,
2530 : // it can falsely identify a DOF at the mid-edge node. This is why
2531 : // we go through FEInterface instead of obj->n_comp() directly.
2532 : const unsigned int nc =
2533 3142916 : FEInterface::n_dofs_at_node(fe_type, &elem, n, add_p_level);
2534 :
2535 : // If this is a non-vertex on a hanging node with extra
2536 : // degrees of freedom, we use the non-vertex dofs (which
2537 : // come in reverse order starting from the end, to
2538 : // simplify p refinement)
2539 3142916 : if (extra_hanging_dofs && nc && !elem.is_vertex(n))
2540 : {
2541 533599 : const int dof_offset = n_comp - nc;
2542 :
2543 : // We should never have fewer dofs than necessary on a
2544 : // node unless we're getting indices on a parent element,
2545 : // and we should never need the indices on such a node
2546 533599 : if (dof_offset < 0)
2547 : {
2548 0 : libmesh_assert(!elem.active());
2549 0 : di.resize(di.size() + nc, DofObject::invalid_id);
2550 : }
2551 : else
2552 1580697 : for (unsigned int i = dof_offset; i != n_comp; ++i)
2553 : {
2554 : const dof_id_type d =
2555 1047098 : obj.dof_number(sys_num, vg, vig, i, n_comp);
2556 344839 : libmesh_assert_not_equal_to (d, DofObject::invalid_id);
2557 1047098 : di.push_back(d);
2558 : }
2559 : }
2560 : // If this is a vertex or an element without extra hanging
2561 : // dofs, our dofs come in forward order coming from the
2562 : // beginning. But we still might not have all those dofs, in cases
2563 : // where a subdomain-restricted variable just had its subdomain
2564 : // expanded.
2565 : else
2566 : {
2567 : const unsigned int good_nc =
2568 2628472 : std::min(static_cast<unsigned int>(n_comp), nc);
2569 5160862 : for (unsigned int i=0; i != good_nc; ++i)
2570 : {
2571 : const dof_id_type d =
2572 2551545 : obj.dof_number(sys_num, vg, vig, i, n_comp);
2573 748003 : libmesh_assert_not_equal_to (d, DofObject::invalid_id);
2574 2551545 : di.push_back(d);
2575 : }
2576 2609317 : for (unsigned int i=good_nc; i != nc; ++i)
2577 0 : di.push_back(DofObject::invalid_id);
2578 : }
2579 3142916 : }
2580 :
2581 : void
2582 45265741 : DofMap::_dof_indices(const Elem & elem,
2583 : int p_level,
2584 : std::vector<dof_id_type> & di,
2585 : const unsigned int vg,
2586 : const unsigned int vig,
2587 : const Node * const * nodes,
2588 : unsigned int n_nodes,
2589 : const unsigned int v
2590 : #ifdef DEBUG
2591 : ,
2592 : std::size_t & tot_size
2593 : #endif
2594 : ) const
2595 : {
2596 45265741 : _dof_indices(elem,
2597 : p_level,
2598 : di,
2599 : vg,
2600 : vig,
2601 : nodes,
2602 : n_nodes,
2603 : v,
2604 : #ifdef DEBUG
2605 : tot_size,
2606 : #endif
2607 : [](const Elem &,
2608 : unsigned int,
2609 : unsigned int,
2610 : std::vector<dof_id_type> & functor_di,
2611 162619038 : const dof_id_type dof) { functor_di.push_back(dof); });
2612 45265741 : }
2613 :
2614 670757 : void DofMap::SCALAR_dof_indices (std::vector<dof_id_type> & di,
2615 : const unsigned int vn,
2616 : #ifdef LIBMESH_ENABLE_AMR
2617 : const bool old_dofs
2618 : #else
2619 : const bool
2620 : #endif
2621 : ) const
2622 : {
2623 : // dof_indices() is a relatively light-weight function; the cost of
2624 : // the logging code itself is roughly on par with the time required
2625 : // to call dof_indices().
2626 : // LOG_SCOPE("SCALAR_dof_indices()", "DofMap");
2627 :
2628 177429 : libmesh_assert(this->variable(vn).type().family == SCALAR);
2629 :
2630 : #ifdef LIBMESH_ENABLE_AMR
2631 : // If we're asking for old dofs then we'd better have some
2632 177429 : if (old_dofs)
2633 66 : libmesh_assert_greater_equal(n_old_dofs(), n_SCALAR_dofs());
2634 :
2635 848186 : dof_id_type my_idx = old_dofs ?
2636 670757 : this->_first_old_scalar_df[vn] : this->_first_scalar_df[vn];
2637 : #else
2638 : dof_id_type my_idx = this->_first_scalar_df[vn];
2639 : #endif
2640 :
2641 177429 : libmesh_assert_not_equal_to(my_idx, DofObject::invalid_id);
2642 :
2643 : // The number of SCALAR dofs comes from the variable order
2644 670757 : const int n_dofs_vn = this->variable(vn).type().order.get_order();
2645 :
2646 670757 : di.resize(n_dofs_vn);
2647 1341514 : for (int i = 0; i != n_dofs_vn; ++i)
2648 848186 : di[i] = my_idx++;
2649 670757 : }
2650 :
2651 :
2652 :
2653 225164 : bool DofMap::semilocal_index (dof_id_type dof_index) const
2654 : {
2655 : // If it's not in the local indices
2656 225164 : if (!this->local_index(dof_index))
2657 : {
2658 : // and if it's not in the ghost indices, then we're not
2659 : // semilocal
2660 188894 : if (!std::binary_search(_send_list.begin(), _send_list.end(), dof_index))
2661 4660 : return false;
2662 : }
2663 :
2664 118350 : return true;
2665 : }
2666 :
2667 :
2668 :
2669 19968 : bool DofMap::all_semilocal_indices (const std::vector<dof_id_type> & dof_indices_in) const
2670 : {
2671 : // We're all semilocal unless we find a counterexample
2672 238160 : for (const auto & di : dof_indices_in)
2673 225164 : if (!this->semilocal_index(di))
2674 2330 : return false;
2675 :
2676 6610 : return true;
2677 : }
2678 :
2679 :
2680 :
2681 : template <typename DofObjectSubclass>
2682 49920 : bool DofMap::is_evaluable(const DofObjectSubclass & obj,
2683 : unsigned int var_num) const
2684 : {
2685 : // Everything is evaluable on a local object
2686 60809 : if (obj.processor_id() == this->processor_id())
2687 11144 : return true;
2688 :
2689 17832 : std::vector<dof_id_type> di;
2690 :
2691 19884 : if (var_num == libMesh::invalid_uint)
2692 1632 : this->dof_indices(&obj, di);
2693 : else
2694 18252 : this->dof_indices(&obj, di, var_num);
2695 :
2696 19884 : return this->all_semilocal_indices(di);
2697 : }
2698 :
2699 :
2700 :
2701 : #ifdef LIBMESH_ENABLE_AMR
2702 :
2703 3625697 : void DofMap::old_dof_indices (const Elem * const elem,
2704 : std::vector<dof_id_type> & di,
2705 : const unsigned int vn) const
2706 : {
2707 : // dof_indices() is a relatively light-weight function; the cost of
2708 : // the logging code itself is roughly on par with the time required
2709 : // to call dof_indices().
2710 : // LOG_SCOPE("old_dof_indices()", "DofMap");
2711 :
2712 1121954 : libmesh_assert(elem);
2713 :
2714 3625697 : const ElemType type = elem->type();
2715 2243369 : const unsigned int sys_num = this->sys_number();
2716 1121954 : const unsigned int n_var_groups = this->n_variable_groups();
2717 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2718 2643530 : const bool is_inf = elem->infinite();
2719 : #endif
2720 :
2721 : // If we have dof indices stored on the elem, and there's no chance
2722 : // that we only have those indices because we were just p refined,
2723 : // then we should have old dof indices too.
2724 1121954 : libmesh_assert(!elem->has_dofs(sys_num) ||
2725 : elem->p_refinement_flag() == Elem::JUST_REFINED ||
2726 : elem->get_old_dof_object());
2727 :
2728 : // Clear the DOF indices vector.
2729 1121954 : di.clear();
2730 :
2731 : // Determine the nodes contributing to element elem
2732 2243908 : std::vector<const Node *> elem_nodes;
2733 : const Node * const * nodes_ptr;
2734 : unsigned int n_nodes;
2735 3625697 : if (elem->type() == TRI3SUBDIVISION)
2736 : {
2737 : // Subdivision surface FE require the 1-ring around elem
2738 0 : const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
2739 0 : MeshTools::Subdivision::find_one_ring(sd_elem, elem_nodes);
2740 0 : nodes_ptr = elem_nodes.data();
2741 0 : n_nodes = cast_int<unsigned int>(elem_nodes.size());
2742 : }
2743 : else
2744 : {
2745 : // All other FE use only the nodes of elem itself
2746 2243369 : nodes_ptr = elem->get_nodes();
2747 3625697 : n_nodes = elem->n_nodes();
2748 : }
2749 :
2750 : // Get the dof numbers
2751 7418009 : for (unsigned int vg=0; vg<n_var_groups; vg++)
2752 : {
2753 1169978 : const VariableGroup & var = this->variable_group(vg);
2754 1169978 : const unsigned int vars_in_group = var.n_variables();
2755 :
2756 7662528 : for (unsigned int vig=0; vig<vars_in_group; vig++)
2757 : {
2758 2391214 : const unsigned int v = var.number(vig);
2759 3870216 : if ((vn == v) || (vn == libMesh::invalid_uint))
2760 : {
2761 3681000 : if (var.type().family == SCALAR &&
2762 160 : (!elem ||
2763 224 : var.active_on_subdomain(elem->subdomain_id())))
2764 : {
2765 : // We asked for this variable, so add it to the vector.
2766 64 : std::vector<dof_id_type> di_new;
2767 224 : this->SCALAR_dof_indices(di_new,v,true);
2768 224 : di.insert( di.end(), di_new.begin(), di_new.end());
2769 : }
2770 : else
2771 3680712 : if (var.active_on_subdomain(elem->subdomain_id()))
2772 : { // Do this for all the variables if one was not specified
2773 : // or just for the specified variable
2774 :
2775 3666906 : FEType fe_type = var.type();
2776 3666906 : const bool add_p_level = fe_type.p_refinement;
2777 :
2778 : // Increase the polynomial order on p refined elements,
2779 : // but make sure you get the right polynomial order for
2780 : // the OLD degrees of freedom
2781 1134664 : int p_adjustment = 0;
2782 4801102 : if (elem->p_refinement_flag() == Elem::JUST_REFINED)
2783 : {
2784 2688 : libmesh_assert_greater (elem->p_level(), 0);
2785 2688 : p_adjustment = -1;
2786 : }
2787 3656799 : else if (elem->p_refinement_flag() == Elem::JUST_COARSENED)
2788 : {
2789 51 : p_adjustment = 1;
2790 : }
2791 3666906 : p_adjustment *= add_p_level;
2792 :
2793 : // Compute the net amount of "extra" order, including Elem::p_level()
2794 3666906 : int extra_order = int(add_p_level*elem->p_level()) + p_adjustment;
2795 :
2796 : const bool extra_hanging_dofs =
2797 3666906 : FEInterface::extra_hanging_dofs(fe_type);
2798 :
2799 : const FEInterface::n_dofs_at_node_ptr ndan =
2800 3666906 : FEInterface::n_dofs_at_node_function(fe_type, elem);
2801 :
2802 : // Get the node-based DOF numbers
2803 30093644 : for (unsigned int n=0; n<n_nodes; n++)
2804 : {
2805 26426738 : const Node * node = nodes_ptr[n];
2806 7774196 : const DofObject & old_dof_obj = node->get_old_dof_object_ref();
2807 :
2808 : // There is a potential problem with h refinement. Imagine a
2809 : // quad9 that has a linear FE on it. Then, on the hanging side,
2810 : // it can falsely identify a DOF at the mid-edge node. This is why
2811 : // we call FEInterface instead of node->n_comp() directly.
2812 : const unsigned int nc =
2813 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2814 19283732 : is_inf ?
2815 29568 : FEInterface::n_dofs_at_node(var.type(), extra_order, elem, n) :
2816 : #endif
2817 34175040 : ndan (type, var.type().order + extra_order, n);
2818 :
2819 26426738 : const int n_comp = old_dof_obj.n_comp_group(sys_num,vg);
2820 :
2821 : // If this is a non-vertex on a hanging node with extra
2822 : // degrees of freedom, we use the non-vertex dofs (which
2823 : // come in reverse order starting from the end, to
2824 : // simplify p refinement)
2825 26426738 : if (extra_hanging_dofs && !elem->is_vertex(n))
2826 : {
2827 2270996 : const int dof_offset = n_comp - nc;
2828 :
2829 : // We should never have fewer dofs than necessary on a
2830 : // node unless we're getting indices on a parent element
2831 : // or a just-coarsened element
2832 2270996 : if (dof_offset < 0)
2833 : {
2834 0 : libmesh_assert(!elem->active() || elem->refinement_flag() ==
2835 : Elem::JUST_COARSENED);
2836 0 : di.resize(di.size() + nc, DofObject::invalid_id);
2837 : }
2838 : else
2839 6155091 : for (int i=n_comp-1; i>=dof_offset; i--)
2840 : {
2841 : const dof_id_type d =
2842 3884095 : old_dof_obj.dof_number(sys_num, vg, vig, i, n_comp);
2843 :
2844 : // On a newly-expanded subdomain, we
2845 : // may have some DoFs that didn't
2846 : // exist in the old system, in which
2847 : // case we can't assert this:
2848 : // libmesh_assert_not_equal_to (d, DofObject::invalid_id);
2849 :
2850 3884095 : di.push_back(d);
2851 : }
2852 : }
2853 : // If this is a vertex or an element without extra hanging
2854 : // dofs, our dofs come in forward order coming from the
2855 : // beginning. But we still might not have all
2856 : // those dofs on the old_dof_obj, in cases
2857 : // where a subdomain-restricted variable just
2858 : // had its subdomain expanded.
2859 : else
2860 : {
2861 : const unsigned int old_nc =
2862 24226827 : std::min(static_cast<unsigned int>(n_comp), nc);
2863 48376768 : for (unsigned int i=0; i != old_nc; ++i)
2864 : {
2865 : const dof_id_type d =
2866 24221026 : old_dof_obj.dof_number(sys_num, vg, vig, i, n_comp);
2867 :
2868 7078297 : libmesh_assert_not_equal_to (d, DofObject::invalid_id);
2869 :
2870 24221026 : di.push_back(d);
2871 : }
2872 24155742 : for (unsigned int i=old_nc; i != nc; ++i)
2873 0 : di.push_back(DofObject::invalid_id);
2874 : }
2875 : }
2876 :
2877 : // If there are any element-based DOF numbers, get them
2878 : const unsigned int nc =
2879 3666906 : FEInterface::n_dofs_per_elem(fe_type, extra_order, elem);
2880 :
2881 3666906 : if (nc != 0)
2882 : {
2883 160653 : const DofObject & old_dof_obj = elem->get_old_dof_object_ref();
2884 :
2885 : const unsigned int n_comp =
2886 160653 : old_dof_obj.n_comp_group(sys_num,vg);
2887 :
2888 502558 : if (old_dof_obj.n_systems() > sys_num &&
2889 : nc <= n_comp)
2890 : {
2891 :
2892 1935936 : for (unsigned int i=0; i<nc; i++)
2893 : {
2894 : const dof_id_type d =
2895 1433378 : old_dof_obj.dof_number(sys_num, vg, vig, i, n_comp);
2896 :
2897 1433378 : di.push_back(d);
2898 : }
2899 : }
2900 : else
2901 : {
2902 : // We should never have fewer dofs than
2903 : // necessary on an element unless we're
2904 : // getting indices on a parent element, a
2905 : // just-coarsened element ... or a
2906 : // subdomain-restricted variable with a
2907 : // just-expanded subdomain
2908 : // libmesh_assert(!elem->active() || fe_type.family == LAGRANGE ||
2909 : // elem->refinement_flag() == Elem::JUST_COARSENED);
2910 0 : di.resize(di.size() + nc, DofObject::invalid_id);
2911 : }
2912 : }
2913 : }
2914 : }
2915 : } // end loop over variables within group
2916 : } // end loop over variable groups
2917 3625697 : }
2918 :
2919 : #endif // LIBMESH_ENABLE_AMR
2920 :
2921 :
2922 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
2923 :
2924 3169993 : void DofMap::find_connected_dofs (std::vector<dof_id_type> & elem_dofs) const
2925 : {
2926 : typedef std::set<dof_id_type> RCSet;
2927 :
2928 : // First insert the DOFS we already depend on into the set.
2929 4042590 : RCSet dof_set (elem_dofs.begin(), elem_dofs.end());
2930 :
2931 872597 : bool done = true;
2932 :
2933 : // Next insert any dofs those might be constrained in terms
2934 : // of. Note that in this case we may not be done: Those may
2935 : // in turn depend on others. So, we need to repeat this process
2936 : // in that case until the system depends only on unconstrained
2937 : // degrees of freedom.
2938 20526857 : for (const auto & dof : elem_dofs)
2939 17356864 : if (this->is_constrained_dof(dof))
2940 : {
2941 : // If the DOF is constrained
2942 : DofConstraints::const_iterator
2943 632840 : pos = _dof_constraints.find(dof);
2944 :
2945 632840 : libmesh_assert (pos != _dof_constraints.end());
2946 :
2947 632840 : const DofConstraintRow & constraint_row = pos->second;
2948 :
2949 : // adaptive p refinement currently gives us lots of empty constraint
2950 : // rows - we should optimize those DoFs away in the future. [RHS]
2951 : //libmesh_assert (!constraint_row.empty());
2952 :
2953 : // Add the DOFs this dof is constrained in terms of.
2954 : // note that these dofs might also be constrained, so
2955 : // we will need to call this function recursively.
2956 6833021 : for (const auto & pr : constraint_row)
2957 3284900 : if (!dof_set.count (pr.first))
2958 : {
2959 1119811 : dof_set.insert (pr.first);
2960 322133 : done = false;
2961 : }
2962 : }
2963 :
2964 :
2965 : // If not done then we need to do more work
2966 : // (obviously :-) )!
2967 3169993 : if (!done)
2968 : {
2969 : // Fill the vector with the contents of the set
2970 125761 : elem_dofs.clear();
2971 311967 : elem_dofs.insert (elem_dofs.end(),
2972 377283 : dof_set.begin(), dof_set.end());
2973 :
2974 :
2975 : // May need to do this recursively. It is possible
2976 : // that we just replaced a constrained DOF with another
2977 : // constrained DOF.
2978 437728 : this->find_connected_dofs (elem_dofs);
2979 :
2980 : } // end if (!done)
2981 3169993 : }
2982 :
2983 : #endif // LIBMESH_ENABLE_CONSTRAINTS
2984 :
2985 :
2986 :
2987 0 : void DofMap::print_info(std::ostream & os) const
2988 : {
2989 0 : os << this->get_info();
2990 0 : }
2991 :
2992 :
2993 :
2994 1711 : std::string DofMap::get_info() const
2995 : {
2996 2715 : std::ostringstream os;
2997 :
2998 : // If we didn't calculate the exact sparsity pattern, the threaded
2999 : // sparsity pattern assembly may have just given us an upper bound
3000 : // on sparsity.
3001 502 : const char * may_equal = " <= ";
3002 :
3003 : // If we calculated the exact sparsity pattern, then we can report
3004 : // exact bandwidth figures:
3005 3441 : for (const auto & mat : _matrices)
3006 1730 : if (mat->need_full_sparsity_pattern())
3007 0 : may_equal = " = ";
3008 :
3009 1711 : dof_id_type max_n_nz = 0, max_n_oz = 0;
3010 1711 : long double avg_n_nz = 0, avg_n_oz = 0;
3011 :
3012 1711 : if (_sp)
3013 : {
3014 1490388 : for (const auto & val : _sp->get_n_nz())
3015 : {
3016 1488866 : max_n_nz = std::max(max_n_nz, val);
3017 1488866 : avg_n_nz += val;
3018 : }
3019 :
3020 1522 : std::size_t n_nz_size = _sp->get_n_nz().size();
3021 :
3022 1522 : this->comm().max(max_n_nz);
3023 1522 : this->comm().sum(avg_n_nz);
3024 1522 : this->comm().sum(n_nz_size);
3025 :
3026 1522 : avg_n_nz /= std::max(n_nz_size,std::size_t(1));
3027 :
3028 1490388 : for (const auto & val : _sp->get_n_oz())
3029 : {
3030 1488866 : max_n_oz = std::max(max_n_oz, val);
3031 1488866 : avg_n_oz += val;
3032 : }
3033 :
3034 1522 : std::size_t n_oz_size = _sp->get_n_oz().size();
3035 :
3036 1522 : this->comm().max(max_n_oz);
3037 1522 : this->comm().sum(avg_n_oz);
3038 1522 : this->comm().sum(n_oz_size);
3039 :
3040 1522 : avg_n_oz /= std::max(n_oz_size,std::size_t(1));
3041 : }
3042 :
3043 : os << " DofMap Sparsity\n Average On-Processor Bandwidth"
3044 2213 : << may_equal << avg_n_nz << '\n';
3045 :
3046 : os << " Average Off-Processor Bandwidth"
3047 2213 : << may_equal << avg_n_oz << '\n';
3048 :
3049 : os << " Maximum On-Processor Bandwidth"
3050 2213 : << may_equal << max_n_nz << '\n';
3051 :
3052 : os << " Maximum Off-Processor Bandwidth"
3053 1711 : << may_equal << max_n_oz << std::endl;
3054 :
3055 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
3056 :
3057 1711 : std::size_t n_constraints = 0, max_constraint_length = 0,
3058 1711 : n_rhss = 0;
3059 1711 : long double avg_constraint_length = 0.;
3060 :
3061 207329 : for (const auto & [constrained_dof, row] : _dof_constraints)
3062 : {
3063 : // Only count local constraints, then sum later
3064 205618 : if (!this->local_index(constrained_dof))
3065 45366 : continue;
3066 :
3067 160252 : std::size_t rowsize = row.size();
3068 :
3069 160252 : max_constraint_length = std::max(max_constraint_length,
3070 42776 : rowsize);
3071 160252 : avg_constraint_length += rowsize;
3072 160252 : n_constraints++;
3073 :
3074 42776 : if (_primal_constraint_values.count(constrained_dof))
3075 14679 : n_rhss++;
3076 : }
3077 :
3078 1711 : this->comm().sum(n_constraints);
3079 1711 : this->comm().sum(n_rhss);
3080 1711 : this->comm().sum(avg_constraint_length);
3081 1711 : this->comm().max(max_constraint_length);
3082 :
3083 1209 : os << " DofMap Constraints\n Number of DoF Constraints = "
3084 1711 : << n_constraints;
3085 1711 : if (n_rhss)
3086 : os << '\n'
3087 499 : << " Number of Heterogenous Constraints= " << n_rhss;
3088 1711 : if (n_constraints)
3089 : {
3090 975 : avg_constraint_length /= n_constraints;
3091 :
3092 : os << '\n'
3093 1265 : << " Average DoF Constraint Length= " << avg_constraint_length;
3094 : }
3095 :
3096 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3097 1004 : std::size_t n_node_constraints = 0, max_node_constraint_length = 0,
3098 1004 : n_node_rhss = 0;
3099 1004 : long double avg_node_constraint_length = 0.;
3100 :
3101 36304 : for (const auto & [node, pr] : _node_constraints)
3102 : {
3103 : // Only count local constraints, then sum later
3104 52950 : if (node->processor_id() != this->processor_id())
3105 6724 : continue;
3106 :
3107 14288 : const NodeConstraintRow & row = pr.first;
3108 28576 : std::size_t rowsize = row.size();
3109 :
3110 28576 : max_node_constraint_length = std::max(max_node_constraint_length,
3111 14288 : rowsize);
3112 28576 : avg_node_constraint_length += rowsize;
3113 28576 : n_node_constraints++;
3114 :
3115 14288 : if (pr.second != Point(0))
3116 0 : n_node_rhss++;
3117 : }
3118 :
3119 1004 : this->comm().sum(n_node_constraints);
3120 1004 : this->comm().sum(n_node_rhss);
3121 1004 : this->comm().sum(avg_node_constraint_length);
3122 1004 : this->comm().max(max_node_constraint_length);
3123 :
3124 1004 : os << "\n Number of Node Constraints = " << n_node_constraints;
3125 1004 : if (n_node_rhss)
3126 : os << '\n'
3127 0 : << " Number of Heterogenous Node Constraints= " << n_node_rhss;
3128 1004 : if (n_node_constraints)
3129 : {
3130 200 : avg_node_constraint_length /= n_node_constraints;
3131 100 : os << "\n Maximum Node Constraint Length= " << max_node_constraint_length
3132 : << '\n'
3133 400 : << " Average Node Constraint Length= " << avg_node_constraint_length;
3134 : }
3135 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
3136 :
3137 502 : os << std::endl;
3138 :
3139 : #endif // LIBMESH_ENABLE_CONSTRAINTS
3140 :
3141 2213 : return os.str();
3142 707 : }
3143 :
3144 54 : void DofMap::create_static_condensation(MeshBase & mesh, System & sys)
3145 : {
3146 72 : _sc = std::make_unique<StaticCondensationDofMap>(mesh, sys, *this);
3147 54 : }
3148 :
3149 38627 : void DofMap::reinit_static_condensation()
3150 : {
3151 38627 : if (_sc)
3152 : {
3153 : // Static condensation reduces the space SNES solves on, which excludes constrained dofs when
3154 : // they are eliminated
3155 120 : _sc->set_parent(_constraint_elimination.get());
3156 360 : _sc->reinit();
3157 : }
3158 38627 : }
3159 :
3160 42 : void DofMap::create_constraint_elimination(MeshBase & mesh, System & sys)
3161 : {
3162 56 : _constraint_elimination = std::make_unique<ConstraintEliminationDofMap>(mesh, sys, *this);
3163 42 : }
3164 :
3165 38627 : void DofMap::reinit_constraint_elimination()
3166 : {
3167 38627 : if (_constraint_elimination)
3168 42 : _constraint_elimination->reinit();
3169 38627 : }
3170 :
3171 38661 : unsigned int DofMap::add_variable(System & sys,
3172 : std::string_view var,
3173 : const FEType & type,
3174 : const std::set<subdomain_id_type> * const active_subdomains)
3175 : {
3176 11110 : parallel_object_only(); // Not strictly needed, but the only safe way to keep in sync
3177 :
3178 11110 : libmesh_assert(this->comm().verify(std::string(var)));
3179 11110 : libmesh_assert(this->comm().verify(type));
3180 11110 : libmesh_assert(this->comm().verify((active_subdomains == nullptr)));
3181 :
3182 11110 : if (active_subdomains)
3183 148 : libmesh_assert(this->comm().verify(active_subdomains->size()));
3184 :
3185 : // Make sure the variable isn't there already
3186 : // or if it is, that it's the type we want
3187 46599 : for (auto v : make_range(this->n_vars()))
3188 7950 : if (this->variable_name(v) == var)
3189 : {
3190 0 : if (this->variable_type(v) == type)
3191 : {
3192 : // Check whether the existing variable's active subdomains also matches
3193 : // the incoming variable's active subdomains. If they don't match, then
3194 : // either it is an error by the user or the user is trying to change the
3195 : // subdomain restriction after the variable has already been added, which
3196 : // is not supported.
3197 12 : const Variable & existing_var = this->variable(v);
3198 :
3199 : // Check whether active_subdomains is not provided/empty and the existing_var is
3200 : // implicitly_active()
3201 12 : bool check1 = (!active_subdomains || active_subdomains->empty()) &&
3202 0 : existing_var.implicitly_active();
3203 :
3204 : // Check if the provided active_subdomains is equal to the existing_var's
3205 : // active_subdomains
3206 : bool check2 =
3207 12 : (active_subdomains && (*active_subdomains == existing_var.active_subdomains()));
3208 :
3209 : // If either of these checks passed, then we already have this variable
3210 12 : if (check1 || check2)
3211 0 : return _variables[v].number();
3212 : }
3213 :
3214 0 : libmesh_error_msg("ERROR: incompatible variable "
3215 : << var << " has already been added for this system!");
3216 : }
3217 :
3218 11110 : libmesh_assert(!sys.is_initialized());
3219 :
3220 38649 : if (this->n_variable_groups())
3221 : {
3222 : // Optimize for VariableGroups here - if the user is adding multiple
3223 : // variables of the same FEType and subdomain restriction, catch
3224 : // that here and add them as members of the same VariableGroup.
3225 : //
3226 : // start by setting this flag to whatever the user has requested
3227 : // and then consider the conditions which should negate it.
3228 1984 : bool should_be_in_vg = this->identify_variable_groups();
3229 :
3230 992 : VariableGroup & vg = _variable_groups.back();
3231 :
3232 : // get a pointer to their subdomain restriction, if any.
3233 : const std::set<subdomain_id_type> * const their_active_subdomains(
3234 3372 : vg.implicitly_active() ? nullptr : &vg.active_subdomains());
3235 :
3236 : // Different types?
3237 2380 : if (vg.type() != type)
3238 342 : should_be_in_vg = false;
3239 :
3240 : // they are restricted, we aren't?
3241 3476 : if (their_active_subdomains &&
3242 346 : (!active_subdomains || (active_subdomains && active_subdomains->empty())))
3243 0 : should_be_in_vg = false;
3244 :
3245 : // they aren't restricted, we are?
3246 3372 : if (!their_active_subdomains && (active_subdomains && !active_subdomains->empty()))
3247 2 : should_be_in_vg = false;
3248 :
3249 3372 : if (their_active_subdomains && active_subdomains)
3250 : // restricted to different sets?
3251 346 : if (*their_active_subdomains != *active_subdomains)
3252 52 : should_be_in_vg = false;
3253 :
3254 : // OK, after all that, append the variable to the vg if none of the conditions
3255 : // were violated
3256 3248 : if (should_be_in_vg)
3257 : {
3258 2108 : const unsigned int vn = this->n_vars();
3259 :
3260 1492 : std::string varstr(var);
3261 :
3262 2108 : _variable_numbers[varstr] = vn;
3263 2108 : vg.append(std::move(varstr));
3264 4216 : _variables.push_back(vg(vg.n_variables() - 1));
3265 2108 : const unsigned int vgn = _variable_groups.size() - 1;
3266 2108 : _variable_group_numbers.push_back(vgn);
3267 616 : _var_to_vg.emplace(vn, vgn);
3268 :
3269 1232 : return vn;
3270 : }
3271 : }
3272 :
3273 : // otherwise, fall back to adding a single variable group
3274 36541 : return this->add_variables(
3275 88635 : sys, std::vector<std::string>(1, std::string(var)), type, active_subdomains);
3276 : }
3277 :
3278 36604 : unsigned int DofMap::add_variables(System & sys,
3279 : const std::vector<std::string> & vars,
3280 : const FEType & type,
3281 : const std::set<subdomain_id_type> * const active_subdomains)
3282 : {
3283 10512 : parallel_object_only(); // Not strictly needed, but the only safe way to keep in sync
3284 :
3285 10512 : libmesh_assert(!sys.is_initialized());
3286 :
3287 10512 : libmesh_assert(this->comm().verify(vars.size()));
3288 10512 : libmesh_assert(this->comm().verify(type));
3289 10512 : libmesh_assert(this->comm().verify((active_subdomains == nullptr)));
3290 :
3291 10512 : if (active_subdomains)
3292 96 : libmesh_assert(this->comm().verify(active_subdomains->size()));
3293 :
3294 : // Make sure the variable isn't there already
3295 : // or if it is, that it's the type we want
3296 773257 : for (auto ovar : vars)
3297 : {
3298 210526 : libmesh_assert(this->comm().verify(ovar));
3299 :
3300 739007 : for (auto v : make_range(this->n_vars()))
3301 1634 : if (this->variable_name(v) == ovar)
3302 : {
3303 0 : if (this->variable_type(v) == type)
3304 0 : return _variables[v].number();
3305 :
3306 0 : libmesh_error_msg("ERROR: incompatible variable "
3307 : << ovar << " has already been added for this system!");
3308 : }
3309 : }
3310 :
3311 36604 : if (this->n_variable_groups())
3312 : {
3313 : // Optimize for VariableGroups here - if the user is adding multiple
3314 : // variables of the same FEType and subdomain restriction, catch
3315 : // that here and add them as members of the same VariableGroup.
3316 : //
3317 : // start by setting this flag to whatever the user has requested
3318 : // and then consider the conditions which should negate it.
3319 752 : bool should_be_in_vg = this->identify_variable_groups();
3320 :
3321 376 : VariableGroup & vg = _variable_groups.back();
3322 :
3323 : // get a pointer to their subdomain restriction, if any.
3324 : const std::set<subdomain_id_type> * const their_active_subdomains(
3325 1264 : vg.implicitly_active() ? nullptr : &vg.active_subdomains());
3326 :
3327 : // Different types?
3328 888 : if (vg.type() != type)
3329 342 : should_be_in_vg = false;
3330 :
3331 : // they are restricted, we aren't?
3332 1316 : if (their_active_subdomains &&
3333 176 : (!active_subdomains || (active_subdomains && active_subdomains->empty())))
3334 0 : should_be_in_vg = false;
3335 :
3336 : // they aren't restricted, we are?
3337 1264 : if (!their_active_subdomains && (active_subdomains && !active_subdomains->empty()))
3338 2 : should_be_in_vg = false;
3339 :
3340 1264 : if (their_active_subdomains && active_subdomains)
3341 : // restricted to different sets?
3342 176 : if (*their_active_subdomains != *active_subdomains)
3343 52 : should_be_in_vg = false;
3344 :
3345 : // If after all that none of the conditions were violated,
3346 : // append the variables to the vg and we're done
3347 1140 : if (should_be_in_vg)
3348 : {
3349 0 : unsigned int vn = this->n_vars();
3350 0 : const unsigned int vgn = _variable_groups.size() - 1;
3351 :
3352 0 : for (auto ovar : vars)
3353 : {
3354 0 : vn = this->n_vars();
3355 :
3356 0 : vg.append(ovar);
3357 :
3358 0 : _variables.push_back(vg(vg.n_variables() - 1));
3359 0 : _variable_numbers[ovar] = vn;
3360 0 : _variable_group_numbers.push_back(vgn);
3361 0 : _var_to_vg.emplace(vn, vgn);
3362 : }
3363 0 : return vn;
3364 : }
3365 : }
3366 :
3367 10512 : const unsigned int curr_n_vars = this->n_vars();
3368 :
3369 36604 : const unsigned int next_first_component = this->n_components(sys.get_mesh());
3370 :
3371 : // We weren't able to add to an existing variable group, so
3372 : // add a new variable group to the list
3373 36604 : _variable_groups.push_back(
3374 : (active_subdomains == nullptr)
3375 73208 : ? VariableGroup(&sys, vars, curr_n_vars, next_first_component, type)
3376 : : VariableGroup(&sys, vars, curr_n_vars, next_first_component, type, *active_subdomains));
3377 :
3378 10512 : const VariableGroup & vg(_variable_groups.back());
3379 36604 : const unsigned int vgn = _variable_groups.size() - 1;
3380 :
3381 : // Add each component of the group individually
3382 773257 : for (auto v : make_range(vars.size()))
3383 : {
3384 736653 : const unsigned int vn = curr_n_vars + v;
3385 1262780 : _variables.push_back(vg(v));
3386 947179 : _variable_numbers[vars[v]] = vn;
3387 736653 : _variable_group_numbers.push_back(vgn);
3388 210526 : _var_to_vg.emplace(vn, vgn);
3389 : }
3390 :
3391 10512 : libmesh_assert_equal_to((curr_n_vars + vars.size()), this->n_vars());
3392 :
3393 : // BSK - Defer this now to System::init_data() so we can detect
3394 : // VariableGroups 12/28/2012
3395 : // // Add the variable group to the _dof_map
3396 : // _dof_map->add_variable_group (vg);
3397 :
3398 : // Return the number of the new variable
3399 57628 : return cast_int<unsigned int>(curr_n_vars + vars.size() - 1);
3400 : }
3401 :
3402 56 : unsigned int DofMap::add_variable_array (System & sys,
3403 : const std::vector<std::string> & vars,
3404 : const FEType & type,
3405 : const std::set<subdomain_id_type> * const active_subdomains)
3406 : {
3407 32 : const unsigned int count = cast_int<unsigned int>(vars.size());
3408 56 : const unsigned int last_var = this->add_variables(sys, vars, type, active_subdomains);
3409 56 : const unsigned int first_var = last_var + 1 - count;
3410 56 : _array_variables.push_back({first_var, first_var + count});
3411 56 : return last_var;
3412 : }
3413 :
3414 45 : void DofMap::get_all_variable_numbers(std::vector<unsigned int> & all_variable_numbers) const
3415 : {
3416 45 : all_variable_numbers.resize(n_vars());
3417 :
3418 14 : unsigned int count = 0;
3419 114 : for (auto vn : _variable_numbers)
3420 91 : all_variable_numbers[count++] = vn.second;
3421 45 : }
3422 :
3423 : template LIBMESH_EXPORT bool DofMap::is_evaluable<Elem>(const Elem &, unsigned int) const;
3424 : template LIBMESH_EXPORT bool DofMap::is_evaluable<Node>(const Node &, unsigned int) const;
3425 :
3426 : } // namespace libMesh
|
