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