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