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