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