Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : #ifndef LIBMESH_DOF_MAP_H
21 : #define LIBMESH_DOF_MAP_H
22 :
23 : // Local Includes
24 : #include "libmesh/libmesh_common.h"
25 : #include "libmesh/reference_counted_object.h"
26 : #include "libmesh/libmesh.h" // libMesh::invalid_uint
27 : #include "libmesh/variable.h"
28 : #include "libmesh/threads.h"
29 : #include "libmesh/threads_allocators.h"
30 : #include "libmesh/elem_range.h"
31 : #include "libmesh/ghosting_functor.h"
32 : #include "libmesh/sparsity_pattern.h"
33 : #include "libmesh/parallel_object.h"
34 : #include "libmesh/point.h"
35 : #include "libmesh/utility.h"
36 : #include "libmesh/elem.h"
37 : #include "libmesh/fe_interface.h"
38 : #include "libmesh/libmesh_logging.h"
39 : #include "libmesh/enum_elem_type.h"
40 : #include "libmesh/mesh_subdivision_support.h"
41 : #include "libmesh/dof_map_base.h"
42 :
43 : // TIMPI includes
44 : #include "timpi/parallel_implementation.h"
45 : #include "timpi/parallel_sync.h"
46 :
47 : // C++ Includes
48 : #include <algorithm>
49 : #include <cstddef>
50 : #include <iterator>
51 : #include <map>
52 : #include <string>
53 : #include <vector>
54 : #include <memory>
55 :
56 : namespace libMesh
57 : {
58 :
59 : // Forward Declarations
60 : class CouplingMatrix;
61 : class DefaultCoupling;
62 : class DirichletBoundary;
63 : class DirichletBoundaries;
64 : class DofMap;
65 : class DofObject;
66 : class FEType;
67 : class MeshBase;
68 : class PeriodicBoundaryBase;
69 : class PeriodicBoundaries;
70 : class System;
71 : class NonlinearImplicitSystem;
72 : class StaticCondensationDofMap;
73 : class ReducedDofMap;
74 : class ConstraintEliminationDofMap;
75 : template <typename T> class DenseVectorBase;
76 : template <typename T> class DenseVector;
77 : template <typename T> class DenseMatrix;
78 : template <typename T> class SparseMatrix;
79 : template <typename T> class NumericVector;
80 : enum Order : int;
81 :
82 :
83 :
84 : // ------------------------------------------------------------
85 : // Do we need constraints for anything?
86 :
87 : #if defined(LIBMESH_ENABLE_AMR) || \
88 : defined(LIBMESH_ENABLE_PERIODIC) || \
89 : defined(LIBMESH_ENABLE_DIRICHLET)
90 : # define LIBMESH_ENABLE_CONSTRAINTS 1
91 : #endif
92 :
93 : // ------------------------------------------------------------
94 : // AMR constraint matrix types
95 :
96 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
97 : /**
98 : * A row of the Dof constraint matrix.
99 : */
100 : typedef std::map<dof_id_type, Real,
101 : std::less<dof_id_type>,
102 : Threads::scalable_allocator<std::pair<const dof_id_type, Real>>> DofConstraintRow;
103 :
104 : /**
105 : * The constraint matrix storage format.
106 : * We're using a class instead of a typedef to allow forward
107 : * declarations and future flexibility. Don't delete this from
108 : * a pointer-to-std-map; the destructor isn't virtual!
109 : */
110 25326 : class DofConstraints : public std::map<dof_id_type,
111 : DofConstraintRow,
112 : std::less<dof_id_type>,
113 : Threads::scalable_allocator<std::pair<const dof_id_type, DofConstraintRow>>>
114 : {
115 : };
116 :
117 : /**
118 : * Storage for DofConstraint right hand sides for a particular
119 : * problem. Each dof id with a non-zero constraint offset
120 : * stores it in such a structure.
121 : */
122 25318 : class DofConstraintValueMap :
123 : public std::map<dof_id_type, Number,
124 : std::less<dof_id_type>,
125 : Threads::scalable_allocator<std::pair<const dof_id_type, Number>>>
126 : {
127 : };
128 :
129 : /**
130 : * Storage for DofConstraint right hand sides for all adjoint
131 : * problems.
132 : */
133 25318 : class AdjointDofConstraintValues :
134 : public std::map<unsigned int, DofConstraintValueMap,
135 : std::less<unsigned int>,
136 : Threads::scalable_allocator
137 : <std::pair<const unsigned int, DofConstraintValueMap>>>
138 : {
139 : };
140 :
141 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
142 : /**
143 : * A row of the Node constraint mapping. Currently this just
144 : * stores the topology of the constrained Nodes, but for forward
145 : * compatibility we also include coefficients, so we could add
146 : * Lagrange-positioned-node constraints later.
147 : */
148 : typedef std::map<const Node *, Real,
149 : std::less<const Node *>,
150 : Threads::scalable_allocator<std::pair<const Node * const, Real>>> NodeConstraintRow;
151 :
152 : /**
153 : * The Node constraint storage format.
154 : * We're using a class instead of a typedef to allow forward
155 : * declarations and future flexibility. Don't delete this from
156 : * a pointer-to-std-map; the destructor isn't virtual!
157 : */
158 10190 : class NodeConstraints : public std::map<const Node *,
159 : std::pair<NodeConstraintRow,Point>,
160 : std::less<const Node *>,
161 : Threads::scalable_allocator<std::pair<const Node * const, std::pair<NodeConstraintRow,Point>>>>
162 : {
163 : };
164 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
165 :
166 : #endif // LIBMESH_ENABLE_CONSTRAINTS
167 :
168 :
169 :
170 : /**
171 : * This class handles the numbering of degrees of freedom on a mesh.
172 : * For systems of equations the class supports a fixed number of variables.
173 : * The degrees of freedom are numbered such that sequential, contiguous blocks
174 : * belong to distinct processors. This is so that the resulting data
175 : * structures will work well with parallel linear algebra packages.
176 : *
177 : * \author Benjamin S. Kirk
178 : * \date 2002-2007
179 : * \brief Manages the degrees of freedom (DOFs) in a simulation.
180 : */
181 : class DofMap : public DofMapBase,
182 : public ReferenceCountedObject<DofMap>
183 : {
184 : public:
185 :
186 : /**
187 : * Constructor. Requires the number of the system for which we
188 : * will be numbering degrees of freedom & the parent object
189 : * we are contained in, which defines our communication space.
190 : */
191 : explicit
192 : DofMap(const unsigned int sys_number,
193 : MeshBase & mesh);
194 :
195 : /**
196 : * Destructor.
197 : */
198 : ~DofMap();
199 :
200 : /**
201 : * Backwards compatibility for prior AugmentSparsityPattern users.
202 : */
203 : class AugmentSparsityPattern : public SparsityPattern::AugmentSparsityPattern
204 : {};
205 :
206 : /**
207 : * Abstract base class to be used to add user-defined parallel
208 : * degree of freedom couplings.
209 : */
210 : class AugmentSendList
211 : {
212 : public:
213 : virtual ~AugmentSendList () = default;
214 :
215 : /**
216 : * User-defined function to augment the send list.
217 : */
218 : virtual void augment_send_list (std::vector<dof_id_type> & send_list) = 0;
219 : };
220 :
221 : /**
222 : * Additional matrices may be attached to this \p DofMap.
223 : * They are initialized to the same sparsity structure as
224 : * the major matrix.
225 : */
226 : void attach_matrix (SparseMatrix<Number> & matrix);
227 :
228 : /**
229 : * Additional matrices may be be temporarily initialized by this \p
230 : * DofMap.
231 : * They are initialized to the same sparsity structure as
232 : * the major matrix.
233 : */
234 : void update_sparsity_pattern(SparseMatrix<Number> & matrix) const;
235 :
236 : /**
237 : * Matrices should not be attached more than once. We can test for
238 : * an already-attached matrix if necessary using \p is_attached
239 : */
240 : bool is_attached (SparseMatrix<Number> & matrix);
241 :
242 : /**
243 : * Distribute dofs on the current mesh. Also builds the send list for
244 : * processor \p proc_id, which defaults to 0 for ease of use in serial
245 : * applications.
246 : * \returns The total number of DOFs for the System, summed across all procs.
247 : */
248 : std::size_t distribute_dofs (MeshBase &);
249 :
250 : /**
251 : * Computes the sparsity pattern for the matrices corresponding to
252 : * \p proc_id and sends that data to Linear Algebra packages for
253 : * preallocation of sparse matrices.
254 : */
255 : void compute_sparsity (const MeshBase &);
256 :
257 : /**
258 : * Returns true iff a sparsity pattern has already been computed.
259 : */
260 : bool computed_sparsity_already () const;
261 :
262 : /**
263 : * Sets the current policy for constructing sparsity patterns: if
264 : * \p use_constraints is true (for robustness), we explicitly
265 : * account for sparsity entries created by constraint matrix pre-
266 : * and post- application. If \p use_constraints is false (for
267 : * speed), we calculate only the sparsity pattern of an
268 : * unconstrained matrix. This is false by default, because in
269 : * nearly all applications our constraints do not increase the
270 : * number of non-zeros required in a sparse matrix.
271 : */
272 : void set_constrained_sparsity_construction(bool use_constraints);
273 :
274 : /**
275 : * Sets _need_full_sparsity_pattern to true regardless of the
276 : * requirements by matrices
277 : */
278 : void full_sparsity_pattern_needed();
279 :
280 : /**
281 : * Sets _need_ghost_constraints to true regardless of the requirements
282 : * by static condensation
283 : */
284 : void ghost_constraints_needed();
285 :
286 : /**
287 : * Returns true iff the current policy when constructing sparsity
288 : * patterns is to explicitly account for sparsity entries created by
289 : * constraint matrix pre- and post- application.
290 : */
291 : bool constrained_sparsity_construction();
292 :
293 : /**
294 : * Clears the sparsity pattern
295 : */
296 : void clear_sparsity();
297 :
298 : /**
299 : * Remove any default ghosting functor(s). User-added ghosting
300 : * functors will be unaffected.
301 : *
302 : * Unless user-added equivalent ghosting functors exist, removing
303 : * the default coupling functor is only safe for explicit solves,
304 : * and removing the default algebraic ghosting functor is only safe
305 : * for codes where no evaluations on neighbor cells (e.g. no jump
306 : * error estimators) are done.
307 : *
308 : * Defaults can be restored manually via add_default_ghosting(), or
309 : * automatically if clear() returns the DofMap to a default state.
310 : */
311 : void remove_default_ghosting();
312 :
313 : /**
314 : * Add the default functor(s) for coupling and algebraic ghosting.
315 : * User-added ghosting functors will be unaffected.
316 : */
317 : void add_default_ghosting();
318 :
319 : /**
320 : * Iterator type for coupling and algebraic ghosting functor ranges.
321 : * This has changed in the past and may change again; code should
322 : * use auto or the type here.
323 : */
324 : typedef std::vector<GhostingFunctor *>::const_iterator GhostingFunctorIterator;
325 :
326 : /**
327 : * Adds a functor which can specify coupling requirements for
328 : * creation of sparse matrices.
329 : * Degree of freedom pairs which match the elements and variables
330 : * returned by these functors will be added to the sparsity pattern,
331 : * and the degrees of freedom which live on other processors will be
332 : * added to the send_list for use on ghosted vectors, and the
333 : * elements which live on other processors will be ghosted on a
334 : * distributed mesh.
335 : *
336 : * GhostingFunctor memory must be managed by the code which calls
337 : * this function; the GhostingFunctor lifetime is expected to extend
338 : * until either the functor is removed or the DofMap is destructed.
339 : *
340 : * When \p to_mesh is true, the \p coupling_functor is also added to
341 : * our associated mesh, to ensure that coupled elements do not get
342 : * lost during mesh distribution. (if coupled elements were
343 : * *already* lost there's no getting them back after the fact,
344 : * sorry)
345 : *
346 : * If \p to_mesh is false, no change to mesh ghosting is made;
347 : * the Mesh must already have ghosting functor(s) specifying a
348 : * superset of \p coupling_functor or this is a horrible bug.
349 : */
350 : void add_coupling_functor(GhostingFunctor & coupling_functor,
351 : bool to_mesh = true);
352 :
353 : /**
354 : * Adds a functor which can specify coupling requirements for
355 : * creation of sparse matrices.
356 : *
357 : * GhostingFunctor memory when using this method is managed by the
358 : * shared_ptr mechanism.
359 : */
360 : void add_coupling_functor(std::shared_ptr<GhostingFunctor> coupling_functor,
361 : bool to_mesh = true)
362 : { _shared_functors[coupling_functor.get()] = coupling_functor;
363 : this->add_coupling_functor(*coupling_functor, to_mesh); }
364 :
365 : /**
366 : * Removes a functor which was previously added to the set of
367 : * coupling functors, from both this DofMap and from the underlying
368 : * mesh.
369 : */
370 : void remove_coupling_functor(GhostingFunctor & coupling_functor);
371 :
372 : /**
373 : * Beginning of range of coupling functors
374 : */
375 318143 : GhostingFunctorIterator coupling_functors_begin() const
376 651856 : { return _coupling_functors.begin(); }
377 :
378 : /**
379 : * End of range of coupling functors
380 : */
381 318143 : GhostingFunctorIterator coupling_functors_end() const
382 651856 : { return _coupling_functors.end(); }
383 :
384 : /**
385 : * Default coupling functor
386 : */
387 0 : DefaultCoupling & default_coupling() { return *_default_coupling; }
388 :
389 : /**
390 : * Adds a functor which can specify algebraic ghosting requirements
391 : * for use with distributed vectors. Degrees of freedom on other
392 : * processors which match the elements and variables returned by
393 : * these functors will be added to the send_list, and the elements
394 : * on other processors will be ghosted on a distributed mesh, so
395 : * that the elements can always be found and the solutions on them
396 : * will always be evaluable.
397 : *
398 : * GhostingFunctor memory must be managed by the code which calls
399 : * this function; the GhostingFunctor lifetime is expected to extend
400 : * until either the functor is removed or the DofMap is destructed.
401 : *
402 : * When \p to_mesh is true, the \p coupling_functor is also added to
403 : * our associated mesh, to ensure that evaluable elements do not get
404 : * lost during mesh distribution. (if evaluable elements were
405 : * *already* lost there's no getting them back after the fact,
406 : * sorry)
407 : *
408 : * If \p to_mesh is false, no change to mesh ghosting is made;
409 : * the Mesh must already have ghosting functor(s) specifying a
410 : * superset of \p evaluable_functor or this is a horrible bug.
411 : */
412 : void add_algebraic_ghosting_functor(GhostingFunctor & evaluable_functor,
413 : bool to_mesh = true);
414 :
415 : /**
416 : * Adds a functor which can specify algebraic ghosting requirements
417 : * for use with distributed vectors.
418 : *
419 : * GhostingFunctor memory when using this method is managed by the
420 : * shared_ptr mechanism.
421 : */
422 : void add_algebraic_ghosting_functor(std::shared_ptr<GhostingFunctor> evaluable_functor,
423 : bool to_mesh = true)
424 : { _shared_functors[evaluable_functor.get()] = evaluable_functor;
425 : this->add_algebraic_ghosting_functor(*evaluable_functor, to_mesh); }
426 :
427 : /**
428 : * Removes a functor which was previously added to the set of
429 : * algebraic ghosting functors, from both this DofMap and from the
430 : * underlying mesh.
431 : */
432 : void remove_algebraic_ghosting_functor(GhostingFunctor & evaluable_functor);
433 :
434 : /**
435 : * Beginning of range of algebraic ghosting functors
436 : */
437 1764 : GhostingFunctorIterator algebraic_ghosting_functors_begin() const
438 12992 : { return _algebraic_ghosting_functors.begin(); }
439 :
440 : /**
441 : * End of range of algebraic ghosting functors
442 : */
443 1764 : GhostingFunctorIterator algebraic_ghosting_functors_end() const
444 12992 : { return _algebraic_ghosting_functors.end(); }
445 :
446 : /**
447 : * Default algebraic ghosting functor
448 : */
449 0 : DefaultCoupling & default_algebraic_ghosting() { return *_default_evaluating; }
450 :
451 : /**
452 : * Attach an object to use to populate the
453 : * sparsity pattern with extra entries.
454 : *
455 : * Care must be taken that when adding entries they are sorted into the Rows
456 : *
457 : * Further, you _must_ modify n_nz and n_oz properly!
458 : *
459 : * This is an advanced function... use at your own peril!
460 : */
461 : void attach_extra_sparsity_object (SparsityPattern::AugmentSparsityPattern & asp)
462 : {
463 : _augment_sparsity_pattern = &asp;
464 : }
465 :
466 : /**
467 : * Attach a function pointer to use as a callback to populate the
468 : * sparsity pattern with extra entries.
469 : *
470 : * Care must be taken that when adding entries they are sorted into the Rows
471 : *
472 : * Further, you _must_ modify n_nz and n_oz properly!
473 : *
474 : * This is an advanced function... use at your own peril!
475 : */
476 : void attach_extra_sparsity_function(void (*func)(SparsityPattern::Graph & sparsity,
477 : std::vector<dof_id_type> & n_nz,
478 : std::vector<dof_id_type> & n_oz,
479 : void *),
480 : void * context = nullptr)
481 : { _extra_sparsity_function = func; _extra_sparsity_context = context; }
482 :
483 : /**
484 : * Attach an object to populate the send_list with extra entries.
485 : * This should only add to the send list, but no checking is done
486 : * to enforce this behavior.
487 : *
488 : * This is an advanced function... use at your own peril!
489 : */
490 : void attach_extra_send_list_object (DofMap::AugmentSendList & asl)
491 : {
492 : _augment_send_list = &asl;
493 : }
494 :
495 : /**
496 : * Attach a function pointer to use as a callback to populate the
497 : * send_list with extra entries.
498 : */
499 : void attach_extra_send_list_function(void (*func)(std::vector<dof_id_type> &, void *),
500 : void * context = nullptr)
501 : { _extra_send_list_function = func; _extra_send_list_context = context; }
502 :
503 : /**
504 : * Takes the \p _send_list vector (which may have duplicate entries)
505 : * and sorts it. The duplicate entries are then removed, resulting in
506 : * a sorted \p _send_list with unique entries. Also calls any user-provided
507 : * methods for adding to the send list.
508 : */
509 : void prepare_send_list ();
510 :
511 : /**
512 : * Clears the \p _send_list vector. This should be done in order to completely
513 : * rebuild the send_list from scratch rather than merely adding to the existing
514 : * send_list.
515 : */
516 21448 : void clear_send_list ()
517 : {
518 21448 : _send_list.clear();
519 21448 : }
520 :
521 : /**
522 : * Clears the \p _send_list vector and then rebuilds it. This may be needed
523 : * in special situations, for example when an algebraic coupling functor cannot
524 : * be added to the \p DofMap until after it is completely setup. Then this method
525 : * can be used to rebuild the send_list once the algebraic coupling functor is
526 : * added. Note that while this will recommunicate constraints with the updated
527 : * send_list, this does assume no new constraints have been added since the previous
528 : * reinit_constraints call.
529 : */
530 : void reinit_send_list (MeshBase & mesh);
531 :
532 :
533 : /**
534 : * \returns A constant reference to the \p _send_list for this processor.
535 : *
536 : * The \p _send_list contains the global indices of all the
537 : * variables in the global solution vector that influence the
538 : * current processor. This information can be used for gathers at
539 : * each solution step to retrieve solution values needed for
540 : * computation.
541 : */
542 250582 : const std::vector<dof_id_type> & get_send_list() const { return _send_list; }
543 :
544 : /**
545 : * \returns A constant reference to the \p _n_nz list for this processor.
546 : *
547 : * The vector contains the bandwidth of the on-processor coupling for each
548 : * row of the global matrix that the current processor owns. This
549 : * information can be used to preallocate space for a parallel sparse matrix.
550 : */
551 : const std::vector<dof_id_type> & get_n_nz() const
552 : {
553 : libmesh_assert(_sp);
554 : return _sp->get_n_nz();
555 : }
556 :
557 : /**
558 : * \returns A constant reference to the \p _n_oz list for this processor.
559 : *
560 : * The vector contains the bandwidth of the off-processor coupling for each
561 : * row of the global matrix that the current processor owns. This
562 : * information can be used to preallocate space for a parallel sparse matrix.
563 : */
564 : const std::vector<dof_id_type> & get_n_oz() const
565 : {
566 : libmesh_assert(_sp);
567 : return _sp->get_n_oz();
568 : }
569 :
570 :
571 : /**
572 : * \returns A constant pointer to the sparsity pattern stored here,
573 : * once that has been computed. Returns null if no sparsity pattern
574 : * has yet been computed.
575 : *
576 : * If need_full_sparsity_pattern is false, the "sparsity pattern"
577 : * may only own n_nz and n_oz lists.
578 : */
579 710 : const SparsityPattern::Build * get_sparsity_pattern() const
580 : {
581 710 : return _sp.get();
582 : }
583 :
584 : /**
585 : * \returns The number of variables in the system
586 : */
587 : unsigned int n_vars() const;
588 :
589 : /**
590 : * \returns The name of variable \p i.
591 : */
592 : const std::string & variable_name(const unsigned int i) const;
593 :
594 : /**
595 : * \returns The total number of scalar components in the system's
596 : * variables. This will equal \p n_vars() in the case of all
597 : * scalar-valued variables. If vector variables are involved, we
598 : * will need to leverage the \p mesh
599 : */
600 : unsigned int n_components(const MeshBase & mesh) const;
601 :
602 : /**
603 : * \returns \p true when \p VariableGroup structures should be
604 : * automatically identified, \p false otherwise.
605 : */
606 : bool identify_variable_groups () const;
607 :
608 : /**
609 : * Toggle automatic \p VariableGroup identification.
610 : */
611 : void identify_variable_groups (const bool);
612 :
613 : /**
614 : * \returns An index, starting from 0 for the first component of the
615 : * first variable, and incrementing for each component of each
616 : * (potentially vector-valued) variable in the system in order.
617 : * For systems with only scalar-valued variables, this will be the
618 : * same as \p var_num
619 : *
620 : * Irony: currently our only non-scalar-valued variable type is
621 : * SCALAR.
622 : */
623 : unsigned int variable_scalar_number (unsigned int var_num,
624 : unsigned int component) const;
625 :
626 : /**
627 : * \returns The finite element type for variable number \p i.
628 : */
629 : const FEType & variable_type (const unsigned int i) const;
630 :
631 : /**
632 : * \returns The finite element type for variable \p var.
633 : */
634 : const FEType & variable_type (std::string_view var) const;
635 :
636 : /**
637 : * \returns The variable number associated with
638 : * the user-specified variable named \p var.
639 : */
640 : unsigned int variable_number (std::string_view var) const;
641 :
642 : /**
643 : * \returns \p true if a variable named \p var exists in this System
644 : */
645 : bool has_variable(std::string_view var) const;
646 :
647 : /**
648 : * Fills \p all_variable_numbers with all the variable numbers for the
649 : * variables that have been added to this system.
650 : */
651 : void get_all_variable_numbers(std::vector<unsigned int> & all_variable_numbers) const;
652 :
653 : /**
654 : * Adds the variable \p var to the list of variables
655 : * for this system. If \p active_subdomains is either \p nullptr
656 : * (the default) or points to an empty set, then it will be assumed that
657 : * \p var has no subdomain restrictions
658 : *
659 : * \returns The index number for the new variable.
660 : */
661 : unsigned int add_variable (System & sys,
662 : std::string_view var,
663 : const FEType & type,
664 : const std::set<subdomain_id_type> * const active_subdomains = nullptr);
665 :
666 : /**
667 : * Adds the variables \p vars to the list of variables
668 : * for this system. If \p active_subdomains is either \p nullptr
669 : * (the default) or points to an empty set, then it will be assumed that
670 : * the \p vars have no subdomain restrictions
671 : *
672 : * \returns The index number for the last of the new variables.
673 : */
674 : unsigned int add_variables (System & sys,
675 : const std::vector<std::string> & vars,
676 : const FEType & type,
677 : const std::set<subdomain_id_type> * const active_subdomains = nullptr);
678 :
679 : /**
680 : * Adds variables \p vars to the list of variables
681 : * for this system. If \p active_subdomains is either \p nullptr
682 : * (the default) or points to an empty set, then it will be assumed that
683 : * the \p vars have no subdomain restrictions. This API will end up
684 : * calling \p this->add_variables(). However, we will additionally store data
685 : * that can be leveraged by the \p DofMap to build degrees of freedom
686 : * containers corresponding to all the variables in this variable array
687 : *
688 : * An 'array variable' is simply a sequence
689 : * of contiguous variable numbers defined by pair where the first member of the pair
690 : * is the first number in the variable sequence and the second member of the pair is
691 : * the number of the last variable in the sequence plus one. Array variables may be
692 : * used in tandem with variable grouping by downstream code to build optimized physics
693 : * kernels since each variable in the array will have the same shape functions.
694 : *
695 : * \returns The index number for the last of the new variables.
696 : */
697 : unsigned int add_variable_array (System & sys,
698 : const std::vector<std::string> & vars,
699 : const FEType & type,
700 : const std::set<subdomain_id_type> * const active_subdomains = nullptr);
701 :
702 : /**
703 : * Specify whether or not we perform an extra (opt-mode enabled) check
704 : * for constraint loops. If a constraint loop is present then
705 : * the system constraints are not valid, so if \p error_on_constraint_loop
706 : * is true we will throw an error in this case.
707 : *
708 : * \note We previously referred to these types of constraints as
709 : * "cyclic" but that has now been deprecated, and these will now
710 : * instead be referred to as "constraint loops" in libMesh.
711 : */
712 : void set_error_on_cyclic_constraint(bool error_on_cyclic_constraint);
713 : void set_error_on_constraint_loop(bool error_on_constraint_loop);
714 :
715 : /**
716 : * \returns The \p VariableGroup description object for group \p g.
717 : */
718 : const VariableGroup & variable_group (const unsigned int c) const;
719 :
720 : const Variable & variable (const unsigned int c) const override;
721 :
722 : /**
723 : * \returns The approximation order for variable \p c.
724 : */
725 : Order variable_order (const unsigned int c) const;
726 :
727 : /**
728 : * \returns The approximation order for \p VariableGroup \p vg.
729 : */
730 : Order variable_group_order (const unsigned int vg) const;
731 :
732 : /**
733 : * \returns The finite element type for \p VariableGroup \p vg.
734 : */
735 : const FEType & variable_group_type (const unsigned int vg) const;
736 :
737 : /**
738 : * \returns The number of variables in the global solution vector. Defaults
739 : * to 1, should be 1 for a scalar equation, 3 for 2D incompressible Navier
740 : * Stokes (u,v,p), etc...
741 : */
742 12236176 : unsigned int n_variable_groups() const
743 24262917 : { return cast_int<unsigned int>(_variable_groups.size()); }
744 :
745 267319 : unsigned int n_variables() const override
746 347113 : { return cast_int<unsigned int>(_variables.size()); }
747 :
748 : /**
749 : * \returns The variable group number that the provided variable number belongs to
750 : */
751 : unsigned int var_group_from_var_number(unsigned int var_num) const;
752 :
753 : /**
754 : * \returns \p true if the variables are capable of being stored in a blocked
755 : * form. Presently, this means that there can only be one variable group,
756 : * and that the group has more than one variable.
757 : */
758 2608 : bool has_blocked_representation() const
759 : {
760 3962 : return ((this->n_variable_groups() == 1) && (this->n_variables() > 1));
761 : }
762 :
763 : /**
764 : * \returns The block size, if the variables are amenable to block storage.
765 : * Otherwise 1.
766 : * This routine was originally designed to enable a blocked storage, but
767 : * it turns out this information is still super useful for solvers even when
768 : * we do not use the blocked storage (e.g., MATMPIBAIJ in PETSc). For example (in PCHMG),
769 : * for a system of PDEs, to construct an efficient multilevel preconditioner, we coarsen
770 : * the matrix of one single PDE instead of the entire huge matrix. In order to
771 : * accomplish this, we need to know many PDEs we have. Another use case,
772 : * the fieldsplit preconditioner can be constructed in place with this info without
773 : * involving any user efforts.
774 : */
775 3962 : unsigned int block_size() const
776 : {
777 2746 : return (this->has_blocked_representation() ? this->n_variables() : 1);
778 : }
779 :
780 : using DofMapBase::n_dofs;
781 : /**
782 : * \returns The total number of degrees of freedom for a particular
783 : * variable \p vn.
784 : */
785 : dof_id_type n_dofs(const unsigned int vn) const
786 : {
787 : dof_id_type n = this->n_local_dofs(vn);
788 : this->comm().sum(n);
789 : return n;
790 : }
791 :
792 : /**
793 : * \returns The number of SCALAR dofs.
794 : */
795 27719 : dof_id_type n_SCALAR_dofs() const { return _n_SCALAR_dofs; }
796 :
797 : using DofMapBase::n_local_dofs;
798 : /**
799 : * \returns The number of degrees of freedom on this processor for a
800 : * particular variable \p vn. This is an O(N) operation on serialized or
801 : * O(N/Nproc) operation on distributed meshes.
802 : */
803 0 : dof_id_type n_local_dofs(const unsigned int vn) const
804 : {
805 : dof_id_type n;
806 0 : this->local_variable_indices(n, _mesh, vn);
807 0 : return n;
808 : }
809 :
810 : /**
811 : * \returns The number of degrees of freedom on each partition for a
812 : * particular variable \p vn.
813 : */
814 36 : std::vector<dof_id_type> n_dofs_per_processor(const unsigned int vn) const
815 : {
816 36 : std::vector<dof_id_type> n_local_dofs(this->n_processors(), 0);
817 36 : this->comm().allgather(this->n_local_dofs(vn), n_local_dofs);
818 36 : return n_local_dofs;
819 : }
820 :
821 : /**
822 : * \returns The processor id that owns the dof index \p dof
823 : */
824 22827 : processor_id_type dof_owner(const dof_id_type dof) const
825 : { std::vector<dof_id_type>::const_iterator ub =
826 22827 : std::upper_bound(_end_df.begin(), _end_df.end(), dof);
827 7471 : libmesh_assert (ub != _end_df.end());
828 22827 : return cast_int<processor_id_type>(ub - _end_df.begin());
829 : }
830 :
831 : void dof_indices (const Elem * const elem,
832 : std::vector<dof_id_type> & di) const;
833 :
834 : /**
835 : * Fills the vector \p di with the global degree of freedom indices
836 : * for the element. For one variable, and potentially for a
837 : * non-default element p refinement level
838 : */
839 : void dof_indices (const Elem * const elem,
840 : std::vector<dof_id_type> & di,
841 : const unsigned int vn,
842 : int p_level = -12345) const override;
843 :
844 : /**
845 : * Fills the vector \p di with the global degree of freedom indices
846 : * for the element. This will aggregate all the degrees of the freedom
847 : * from the variable array that \p vn is a member of, and potentially for a
848 : * non-default element p refinement level
849 : */
850 : void array_dof_indices(const Elem * const elem,
851 : std::vector<dof_id_type> & di,
852 : const unsigned int vn,
853 : int p_level = -12345) const;
854 :
855 : void array_dof_indices(const Node * const node,
856 : std::vector<dof_id_type> & di,
857 : const unsigned int vn) const;
858 :
859 : template <typename DofIndicesFunctor>
860 : void array_dof_indices(const DofIndicesFunctor & functor,
861 : std::vector<dof_id_type> & di,
862 : const unsigned int vn) const;
863 :
864 : /**
865 : * Retrieves degree of freedom indices for a given \p elem and then performs actions for these
866 : * indices defined by the user-provided functors \p scalar_dofs_functor and \p field_dofs_functor.
867 : * This API is useful when a user wants to do more than simply fill a degree of freedom container
868 : * @param elem The element to get degrees of freedom for
869 : * @param di A container for degrees of freedom. It is up to the provided functors how this gets
870 : * filled
871 : * @param vn The variable number to retrieve degrees of freedom for
872 : * @param scalar_dofs_functor The functor that acts on scalar degrees of freedom. This functor has
873 : * the interface:
874 : * void scalar_dofs_functor(const Elem & elem,
875 : * std::vector<dof_id_type> & di,
876 : * const std::vector<dof_id_type> & scalar_dof_indices)
877 : * where \p di is the degree of freedom container described above and
878 : * \p scalar_dof_indices are the scalar dof indices available to
879 : * \p elem
880 : * @param field_dofs_functor The functor that acts on "field" (e.g. non-scalar, non-global)
881 : * degrees of freedom. This functor has
882 : * the interface:
883 : * void field_dofs_functor(const Elem & elem,
884 : * const unsigned int node_num,
885 : * const unsigned int var_num,
886 : * std::vector<dof_id_type> & di,
887 : * const dof_id_type field_dof)
888 : * where \p field_dof represents a field degree of freedom to act on and
889 : * is associated with \p node_num and \p var_num. If the degree of
890 : * freedom is elemental than \p node_num will be \p invalid_uint. \p di
891 : * is again the degree of freedom container provided above
892 : */
893 : template <typename ScalarDofsFunctor, typename FieldDofsFunctor>
894 : void dof_indices(const Elem * const elem,
895 : std::vector<dof_id_type> & di,
896 : const unsigned int vn,
897 : ScalarDofsFunctor scalar_dofs_functor,
898 : FieldDofsFunctor field_dofs_functor,
899 : int p_level = -12345) const;
900 :
901 : /**
902 : * Fills the vector \p di with the global degree of freedom indices
903 : * for the \p node.
904 : */
905 : void dof_indices (const Node * const node,
906 : std::vector<dof_id_type> & di) const;
907 :
908 : /**
909 : * Fills the vector \p di with the global degree of freedom indices
910 : * for the \p node, for one variable \p vn.
911 : */
912 : void dof_indices (const Node * const node,
913 : std::vector<dof_id_type> & di,
914 : const unsigned int vn) const override;
915 :
916 : /**
917 : * Appends to the vector \p di the global degree of freedom indices
918 : * for \p elem.node_ref(n), for one variable \p vn. On hanging
919 : * nodes with both vertex and non-vertex DoFs, only those indices
920 : * which are directly supported on \p elem are included.
921 : */
922 : void dof_indices (const Elem & elem,
923 : unsigned int n,
924 : std::vector<dof_id_type> & di,
925 : const unsigned int vn) const;
926 :
927 : #ifdef LIBMESH_ENABLE_AMR
928 :
929 : /**
930 : * Appends to the vector \p di the old global degree of freedom
931 : * indices for \p elem.node_ref(n), for one variable \p vn. On
932 : * hanging nodes with both vertex and non-vertex DoFs, only those
933 : * indices which are directly supported on \p elem are included.
934 : */
935 : void old_dof_indices (const Elem & elem,
936 : unsigned int n,
937 : std::vector<dof_id_type> & di,
938 : const unsigned int vn) const;
939 :
940 : #endif // LIBMESH_ENABLE_AMR
941 :
942 : /**
943 : * Fills the vector \p di with the global degree of freedom indices
944 : * corresponding to the SCALAR variable vn. If old_dofs=true,
945 : * the old SCALAR dof indices are returned.
946 : *
947 : * \note We do not need to pass in an element since SCALARs are
948 : * global variables.
949 : */
950 : void SCALAR_dof_indices (std::vector<dof_id_type> & di,
951 : const unsigned int vn,
952 : const bool old_dofs=false) const;
953 :
954 : /**
955 : * \returns \p true if degree of freedom index \p dof_index
956 : * is either a local index or in the \p send_list.
957 : *
958 : * \note This is an O(logN) operation for a send_list of size N; we
959 : * don't cache enough information for O(1) right now.
960 : */
961 : bool semilocal_index (dof_id_type dof_index) const;
962 :
963 : /**
964 : * \returns \p true if all degree of freedom indices in \p
965 : * dof_indices are either local indices or in the \p send_list.
966 : *
967 : * \note This is an O(logN) operation for a send_list of size N; we
968 : * don't cache enough information for O(1) right now.
969 : */
970 : bool all_semilocal_indices (const std::vector<dof_id_type> & dof_indices) const;
971 :
972 : /**
973 : * \returns \p true if degree of freedom index \p dof_index
974 : * is a local index.
975 : */
976 12949003 : bool local_index (dof_id_type dof_index) const
977 20374884 : { return (dof_index >= this->first_dof()) && (dof_index < this->end_dof()); }
978 :
979 : /**
980 : * \returns \p true iff our solutions can be locally evaluated on
981 : * \p obj (which should be an Elem or a Node) for variable number \p
982 : * var_num (for all variables, if \p var_num is invalid_uint)
983 : */
984 : template <typename DofObjectSubclass>
985 : bool is_evaluable(const DofObjectSubclass & obj,
986 : unsigned int var_num = libMesh::invalid_uint) const;
987 :
988 : /**
989 : * Allow the implicit_neighbor_dofs flag to be set programmatically.
990 : * This overrides the --implicit_neighbor_dofs commandline option.
991 : * We can use this to set the implicit neighbor dofs option differently
992 : * for different systems, whereas the commandline option is the same
993 : * for all systems.
994 : */
995 : void set_implicit_neighbor_dofs(bool implicit_neighbor_dofs);
996 :
997 : /**
998 : * Set the _verify_dirichlet_bc_consistency flag.
999 : */
1000 : void set_verify_dirichlet_bc_consistency(bool val);
1001 :
1002 : /**
1003 : * Tells other library functions whether or not this problem
1004 : * includes coupling between dofs in neighboring cells, as can
1005 : * currently be specified on the command line or inferred from
1006 : * the use of all discontinuous variables.
1007 : */
1008 : bool use_coupled_neighbor_dofs(const MeshBase & mesh) const;
1009 :
1010 : /**
1011 : * Builds the local element vector \p Ue from the global vector \p Ug,
1012 : * accounting for any constrained degrees of freedom. For an element
1013 : * without constrained degrees of freedom this is the trivial mapping
1014 : * \f$ Ue[i] = Ug[dof_indices[i]] \f$
1015 : *
1016 : * \note The user must ensure that the element vector \p Ue is
1017 : * properly sized when calling this method. This is because there
1018 : * is no \p resize() method in the \p DenseVectorBase<> class.
1019 : */
1020 : void extract_local_vector (const NumericVector<Number> & Ug,
1021 : const std::vector<dof_id_type> & dof_indices,
1022 : DenseVectorBase<Number> & Ue) const;
1023 :
1024 : /**
1025 : * If T == dof_id_type, counts, if T == std::vector<dof_id_type>, fills an
1026 : * array of, those dof indices which belong to the given variable number and
1027 : * live on the current processor.
1028 : */
1029 : template <typename T, std::enable_if_t<std::is_same_v<T, dof_id_type> ||
1030 : std::is_same_v<T, std::vector<dof_id_type>>, int> = 0>
1031 : void local_variable_indices(T & idx,
1032 : const MeshBase & mesh,
1033 : unsigned int var_num) const;
1034 :
1035 : /**
1036 : * If T == dof_id_type, counts, if T == std::vector<dof_id_type>, fills an
1037 : * array of, those dof indices which belong to the given variable number and
1038 : * live on the current processor.
1039 : */
1040 : template <typename T,
1041 : std::enable_if_t<std::is_same_v<T, dof_id_type> ||
1042 : std::is_same_v<T, std::vector<dof_id_type>>,
1043 : int> = 0>
1044 : void local_variable_indices(T & idx, unsigned int var_num) const
1045 : { this->local_variable_indices(idx, this->_mesh, var_num); }
1046 :
1047 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
1048 :
1049 : //--------------------------------------------------------------------
1050 : // Constraint-specific methods
1051 : /**
1052 : * \returns The total number of constrained degrees of freedom
1053 : * in the problem.
1054 : */
1055 : dof_id_type n_constrained_dofs() const;
1056 :
1057 : /**
1058 : * \returns The number of constrained degrees of freedom
1059 : * on this processor.
1060 : */
1061 : dof_id_type n_local_constrained_dofs() const;
1062 :
1063 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1064 : /**
1065 : * \returns The total number of constrained Nodes
1066 : * in the mesh.
1067 : */
1068 : dof_id_type n_constrained_nodes() const
1069 : { return cast_int<dof_id_type>(_node_constraints.size()); }
1070 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
1071 :
1072 : /**
1073 : * Rebuilds the raw degree of freedom and DofObject constraints,
1074 : * based on attached DirichletBoundary objects and on non-conforming
1075 : * interface in adapted meshes.
1076 : *
1077 : * A time is specified for use in building time-dependent Dirichlet
1078 : * constraints.
1079 : */
1080 : void create_dof_constraints (const MeshBase &, Real time=0);
1081 :
1082 : /**
1083 : * Gathers constraint equation dependencies from other processors
1084 : */
1085 : void allgather_recursive_constraints (MeshBase &);
1086 :
1087 : /**
1088 : * Sends constraint equations to constraining processors
1089 : */
1090 : void scatter_constraints (MeshBase &);
1091 :
1092 : /**
1093 : * Helper function for querying about constraint equations on other
1094 : * processors. If any id in \p requested_dof_ids is constrained on
1095 : * another processor, its constraint will be added on this processor
1096 : * as well. If \p look_for_constrainees is true, then constraints
1097 : * will also be returned if the id appears as a constraining value
1098 : * not just if it appears as a constrained value.
1099 : *
1100 : * This function operates recursively: if the constraint for a
1101 : * constrained dof is newly added locally, then any other dofs which
1102 : * constrain it are queried to see if they are in turn constrained,
1103 : * and so on.
1104 : */
1105 : void gather_constraints (MeshBase & mesh,
1106 : std::set<dof_id_type> & unexpanded_dofs,
1107 : bool look_for_constrainees);
1108 :
1109 : /**
1110 : * Postprocesses any constrained degrees of freedom
1111 : * to be constrained only in terms of unconstrained dofs, then adds
1112 : * unconstrained dofs to the send_list and prepares that for use.
1113 : * This should be run after both system (create_dof_constraints) and
1114 : * user constraints have all been added.
1115 : */
1116 : void process_constraints (MeshBase &);
1117 :
1118 : /**
1119 : * Throw an error if we detect any constraint loops, i.e.
1120 : * A -> B -> C -> A
1121 : * that is, "dof A is constrained in terms of dof B which is
1122 : * constrained in terms of dof C which is constrained in terms of
1123 : * dof A", since these are not supported by libMesh and give
1124 : * erroneous results if they are present.
1125 : *
1126 : * \note The original "cyclic constraint" terminology was
1127 : * unfortunate since the word cyclic is used by some software to
1128 : * indicate an actual type of rotational/angular constraint and not
1129 : * (as here) a cyclic graph. The former nomenclature will eventually
1130 : * be deprecated in favor of "constraint loop".
1131 : */
1132 : void check_for_cyclic_constraints();
1133 : void check_for_constraint_loops();
1134 :
1135 : /**
1136 : * Adds a copy of the user-defined row to the constraint matrix, using
1137 : * an inhomogeneous right-hand-side for the constraint equation.
1138 : */
1139 : void add_constraint_row (const dof_id_type dof_number,
1140 : const DofConstraintRow & constraint_row,
1141 : const Number constraint_rhs,
1142 : const bool forbid_constraint_overwrite);
1143 :
1144 : /**
1145 : * Adds a copy of the user-defined row to the constraint matrix,
1146 : * using an inhomogeneous right-hand-side for the adjoint constraint
1147 : * equation.
1148 : *
1149 : * \p forbid_constraint_overwrite here only tests for overwriting
1150 : * the rhs. This method should only be used when an equivalent
1151 : * constraint (with a potentially different rhs) already exists for
1152 : * the primal problem.
1153 : */
1154 : void add_adjoint_constraint_row (const unsigned int qoi_index,
1155 : const dof_id_type dof_number,
1156 : const DofConstraintRow & constraint_row,
1157 : const Number constraint_rhs,
1158 : const bool forbid_constraint_overwrite);
1159 :
1160 : /**
1161 : * Adds a copy of the user-defined row to the constraint matrix, using
1162 : * a homogeneous right-hand-side for the constraint equation.
1163 : * By default, produces an error if the DOF was already constrained.
1164 : */
1165 17640 : void add_constraint_row (const dof_id_type dof_number,
1166 : const DofConstraintRow & constraint_row,
1167 : const bool forbid_constraint_overwrite = true)
1168 76046 : { add_constraint_row(dof_number, constraint_row, 0., forbid_constraint_overwrite); }
1169 :
1170 : /**
1171 : * \returns An iterator pointing to the first DoF constraint row.
1172 : */
1173 : DofConstraints::const_iterator constraint_rows_begin() const
1174 : { return _dof_constraints.begin(); }
1175 :
1176 : /**
1177 : * \returns An iterator pointing just past the last DoF constraint row.
1178 : */
1179 : DofConstraints::const_iterator constraint_rows_end() const
1180 : { return _dof_constraints.end(); }
1181 :
1182 : /**
1183 : * Provide a const accessor to the DofConstraints map. This allows the user
1184 : * to quickly search the data structure rather than just iterating over it.
1185 : */
1186 289594 : const DofConstraints & get_dof_constraints() const { return _dof_constraints; }
1187 :
1188 : void stash_dof_constraints()
1189 : {
1190 : libmesh_assert(_stashed_dof_constraints.empty());
1191 : _dof_constraints.swap(_stashed_dof_constraints);
1192 : }
1193 :
1194 : void unstash_dof_constraints()
1195 : {
1196 : libmesh_assert(_dof_constraints.empty());
1197 : _dof_constraints.swap(_stashed_dof_constraints);
1198 : }
1199 :
1200 : /**
1201 : * Similar to the stash/unstash_dof_constraints() API, but swaps
1202 : * _dof_constraints and _stashed_dof_constraints without asserting
1203 : * that the source or destination is empty first.
1204 : *
1205 : * \note There is an implicit assumption that swapping between sets
1206 : * of Constraints does not change the sparsity pattern or expand the
1207 : * send_list, since the only thing changed is the DofConstraints
1208 : * themselves. This is intended to work for swapping between
1209 : * DofConstraints A and B, where A is used to define the send_list,
1210 : * and B is a subset of A.
1211 : */
1212 : void swap_dof_constraints()
1213 : {
1214 : _dof_constraints.swap(_stashed_dof_constraints);
1215 : }
1216 :
1217 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1218 : /**
1219 : * \returns An iterator pointing to the first Node constraint row.
1220 : */
1221 : NodeConstraints::const_iterator node_constraint_rows_begin() const
1222 : { return _node_constraints.begin(); }
1223 :
1224 : /**
1225 : * \returns An iterator pointing just past the last Node constraint row.
1226 : */
1227 : NodeConstraints::const_iterator node_constraint_rows_end() const
1228 : { return _node_constraints.end(); }
1229 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
1230 :
1231 : /**
1232 : * \returns \p true if the degree of freedom dof is constrained,
1233 : * \p false otherwise.
1234 : */
1235 : bool is_constrained_dof (const dof_id_type dof) const;
1236 :
1237 : /**
1238 : * \returns \p true if the system has any heterogeneous constraints for
1239 : * adjoint solution \p qoi_num, \p false otherwise.
1240 : */
1241 : bool has_heterogeneous_adjoint_constraints (const unsigned int qoi_num) const;
1242 :
1243 : /**
1244 : * Backwards compatibility with misspelling.
1245 : */
1246 3072 : bool has_heterogenous_adjoint_constraints (const unsigned int qoi_num) const
1247 : {
1248 6186 : return this->has_heterogeneous_adjoint_constraints (qoi_num);
1249 : }
1250 :
1251 : /**
1252 : * \returns The heterogeneous constraint value if the degree of
1253 : * freedom \p dof has a heterogeneous constraint for adjoint solution
1254 : * \p qoi_num, zero otherwise.
1255 : */
1256 : Number has_heterogeneous_adjoint_constraint (const unsigned int qoi_num,
1257 : const dof_id_type dof) const;
1258 :
1259 : /**
1260 : * Backwards compatibility with misspelling.
1261 : */
1262 1093446 : Number has_heterogenous_adjoint_constraint (const unsigned int qoi_num,
1263 : const dof_id_type dof) const
1264 : {
1265 1244328 : return this->has_heterogeneous_adjoint_constraint (qoi_num, dof);
1266 : }
1267 :
1268 : /**
1269 : * \returns A reference to the set of right-hand-side values in
1270 : * primal constraint equations
1271 : */
1272 : DofConstraintValueMap & get_primal_constraint_values();
1273 :
1274 : /**
1275 : * \returns \p true if the Node is constrained,
1276 : * false otherwise.
1277 : */
1278 : bool is_constrained_node (const Node * node) const;
1279 :
1280 : /**
1281 : * Prints (from processor 0) all DoF and Node constraints. If \p
1282 : * print_nonlocal is true, then each constraint is printed once for
1283 : * each processor that knows about it, which may be useful for \p
1284 : * DistributedMesh debugging.
1285 : */
1286 : void print_dof_constraints(std::ostream & os=libMesh::out,
1287 : bool print_nonlocal=false) const;
1288 :
1289 : /**
1290 : * Gets a string reporting all DoF and Node constraints local to
1291 : * this processor. If \p print_nonlocal is true, then nonlocal
1292 : * constraints which are locally known are included.
1293 : */
1294 : std::string get_local_constraints(bool print_nonlocal=false) const;
1295 :
1296 :
1297 : /**
1298 : * Tests the constrained degrees of freedom on the numeric vector \p v, which
1299 : * represents a solution defined on the mesh, returning a pair whose first
1300 : * entry is the maximum absolute error on a constrained DoF and whose second
1301 : * entry is the maximum relative error. Useful for debugging purposes.
1302 : *
1303 : * If \p v == nullptr, the system solution vector is tested.
1304 : */
1305 : std::pair<Real, Real> max_constraint_error(const System & system,
1306 : NumericVector<Number> * v = nullptr) const;
1307 :
1308 : #endif // LIBMESH_ENABLE_CONSTRAINTS
1309 :
1310 : //--------------------------------------------------------------------
1311 : // Constraint-specific methods
1312 : // Some of these methods are enabled (but inlined away to nothing)
1313 : // when constraints are disabled at configure-time. This is to
1314 : // increase API compatibility of user code with different library
1315 : // builds.
1316 :
1317 : /**
1318 : * Constrains the element matrix. This method requires the
1319 : * element matrix to be square, in which case the elem_dofs
1320 : * correspond to the global DOF indices of both the rows and
1321 : * columns of the element matrix. For this case the rows
1322 : * and columns of the matrix necessarily correspond to variables
1323 : * of the same approximation order.
1324 : *
1325 : * If \p asymmetric_constraint_rows is set to true (as it is by
1326 : * default), constraint row equations will be reinforced in a way
1327 : * which breaks matrix symmetry but makes inexact linear solver
1328 : * solutions more likely to satisfy hanging node constraints.
1329 : */
1330 : void constrain_element_matrix (DenseMatrix<Number> & matrix,
1331 : std::vector<dof_id_type> & elem_dofs,
1332 : bool asymmetric_constraint_rows = true) const;
1333 :
1334 : /**
1335 : * Constrains the element matrix. This method allows the
1336 : * element matrix to be non-square, in which case the row_dofs
1337 : * and col_dofs may be of different size and correspond to
1338 : * variables approximated in different spaces.
1339 : */
1340 : void constrain_element_matrix (DenseMatrix<Number> & matrix,
1341 : std::vector<dof_id_type> & row_dofs,
1342 : std::vector<dof_id_type> & col_dofs,
1343 : bool asymmetric_constraint_rows = true) const;
1344 :
1345 : /**
1346 : * Constrains the element vector.
1347 : */
1348 : void constrain_element_vector (DenseVector<Number> & rhs,
1349 : std::vector<dof_id_type> & dofs,
1350 : bool asymmetric_constraint_rows = true) const;
1351 :
1352 : /**
1353 : * Constrains the element matrix and vector. This method requires
1354 : * the element matrix to be square, in which case the elem_dofs
1355 : * correspond to the global DOF indices of both the rows and
1356 : * columns of the element matrix. For this case the rows
1357 : * and columns of the matrix necessarily correspond to variables
1358 : * of the same approximation order.
1359 : */
1360 : void constrain_element_matrix_and_vector (DenseMatrix<Number> & matrix,
1361 : DenseVector<Number> & rhs,
1362 : std::vector<dof_id_type> & elem_dofs,
1363 : bool asymmetric_constraint_rows = true) const;
1364 :
1365 : /**
1366 : * Constrains the element matrix and vector. This method requires
1367 : * the element matrix to be square, in which case the elem_dofs
1368 : * correspond to the global DOF indices of both the rows and
1369 : * columns of the element matrix. For this case the rows
1370 : * and columns of the matrix necessarily correspond to variables
1371 : * of the same approximation order.
1372 : *
1373 : * The heterogeneous version of this method creates linear systems in
1374 : * which heterogeneously constrained degrees of freedom will solve to
1375 : * their correct offset values, as would be appropriate for finding
1376 : * a solution to a linear problem in a single algebraic solve. The
1377 : * non-heterogeneous version of this method creates linear systems in
1378 : * which even heterogeneously constrained degrees of freedom are
1379 : * solved without offset values taken into account, as would be
1380 : * appropriate for finding linearized updates to a solution in which
1381 : * heterogeneous constraints are already satisfied.
1382 : *
1383 : * By default, the constraints for the primal solution of this
1384 : * system are used. If a non-negative \p qoi_index is passed in,
1385 : * then the constraints for the corresponding adjoint solution are
1386 : * used instead.
1387 : */
1388 : void heterogeneously_constrain_element_matrix_and_vector (DenseMatrix<Number> & matrix,
1389 : DenseVector<Number> & rhs,
1390 : std::vector<dof_id_type> & elem_dofs,
1391 : bool asymmetric_constraint_rows = true,
1392 : int qoi_index = -1) const;
1393 :
1394 : /*
1395 : * Backwards compatibility with misspelling.
1396 : */
1397 0 : void heterogenously_constrain_element_matrix_and_vector (DenseMatrix<Number> & matrix,
1398 : DenseVector<Number> & rhs,
1399 : std::vector<dof_id_type> & elem_dofs,
1400 : bool asymmetric_constraint_rows = true,
1401 : int qoi_index = -1) const
1402 : {
1403 : return this->heterogeneously_constrain_element_matrix_and_vector
1404 4615 : (matrix, rhs, elem_dofs, asymmetric_constraint_rows, qoi_index);
1405 : }
1406 :
1407 : /**
1408 : * Constrains the element vector. This method requires
1409 : * the element matrix to be square and not-yet-constrained, in which
1410 : * case the elem_dofs correspond to the global DOF indices of both
1411 : * the rows and columns of the element matrix.
1412 : *
1413 : * The heterogeneous version of this method creates linear systems in
1414 : * which heterogeneously constrained degrees of freedom will solve to
1415 : * their correct offset values, as would be appropriate for finding
1416 : * a solution to a linear problem in a single algebraic solve. The
1417 : * non-heterogeneous version of this method creates linear systems in
1418 : * which even heterogeneously constrained degrees of freedom are
1419 : * solved without offset values taken into account, as would be
1420 : * appropriate for finding linearized updates to a solution in which
1421 : * heterogeneous constraints are already satisfied.
1422 : *
1423 : * Note the sign difference from the nonlinear heterogeneous constraint
1424 : * method: Solving u:=K\f has the opposite sign convention from
1425 : * u:=u_in-J\r, and we apply heterogeneous constraints accordingly.
1426 : *
1427 : * By default, the constraints for the primal solution of this
1428 : * system are used. If a non-negative \p qoi_index is passed in,
1429 : * then the constraints for the corresponding adjoint solution are
1430 : * used instead.
1431 : */
1432 : void heterogeneously_constrain_element_vector (const DenseMatrix<Number> & matrix,
1433 : DenseVector<Number> & rhs,
1434 : std::vector<dof_id_type> & elem_dofs,
1435 : bool asymmetric_constraint_rows = true,
1436 : int qoi_index = -1) const;
1437 :
1438 : /*
1439 : * Backwards compatibility with misspelling.
1440 : */
1441 140 : void heterogenously_constrain_element_vector (const DenseMatrix<Number> & matrix,
1442 : DenseVector<Number> & rhs,
1443 : std::vector<dof_id_type> & elem_dofs,
1444 : bool asymmetric_constraint_rows = true,
1445 : int qoi_index = -1) const
1446 : {
1447 : return this->heterogeneously_constrain_element_vector
1448 560 : (matrix, rhs, elem_dofs, asymmetric_constraint_rows, qoi_index);
1449 : }
1450 :
1451 : /**
1452 : * Constrains the element Jacobian and residual. The element
1453 : * Jacobian is square, and the elem_dofs should correspond to the
1454 : * global DOF indices of both the rows and columns of the element
1455 : * matrix.
1456 : *
1457 : * The residual-constraining version of this method creates linear
1458 : * systems in which heterogeneously constrained degrees of freedom
1459 : * create non-zero residual terms when not at their correct offset
1460 : * values, as would be appropriate for finding a solution to a
1461 : * nonlinear problem in a quasi-Newton solve.
1462 : *
1463 : * Note the sign difference from the linear heterogeneous constraint
1464 : * method: Solving u:=u_in-J\r has the opposite sign convention from
1465 : * u:=K\f, and we apply heterogeneous constraints accordingly.
1466 : *
1467 : * The \p solution vector passed in should be a serialized or
1468 : * ghosted primal solution
1469 : */
1470 : void heterogeneously_constrain_element_jacobian_and_residual (DenseMatrix<Number> & matrix,
1471 : DenseVector<Number> & rhs,
1472 : std::vector<dof_id_type> & elem_dofs,
1473 : NumericVector<Number> & solution_local) const;
1474 :
1475 : /**
1476 : * Constrains the element residual. The element Jacobian is square,
1477 : * and the elem_dofs should correspond to the global DOF indices of
1478 : * both the rows and columns of the element matrix.
1479 : *
1480 : * The residual-constraining version of this method creates linear
1481 : * systems in which heterogeneously constrained degrees of freedom
1482 : * create non-zero residual terms when not at their correct offset
1483 : * values, as would be appropriate for finding a solution to a
1484 : * nonlinear problem in a quasi-Newton solve.
1485 : *
1486 : * The \p solution vector passed in should be a serialized or
1487 : * ghosted primal solution
1488 : */
1489 : void heterogeneously_constrain_element_residual (DenseVector<Number> & rhs,
1490 : std::vector<dof_id_type> & elem_dofs,
1491 : NumericVector<Number> & solution_local) const;
1492 :
1493 :
1494 : /**
1495 : * Constrains the element residual. The element Jacobian is square,
1496 : * and the elem_dofs should correspond to the global DOF indices of
1497 : * both the rows and columns of the element matrix, and the dof
1498 : * constraint should not include any heterogeneous terms.
1499 : *
1500 : * The residual-constraining version of this method creates linear
1501 : * systems in which heterogeneously constrained degrees of freedom
1502 : * create non-zero residual terms when not at their correct offset
1503 : * values, as would be appropriate for finding a solution to a
1504 : * nonlinear problem in a quasi-Newton solve.
1505 : *
1506 : * The \p solution vector passed in should be a serialized or
1507 : * ghosted primal solution
1508 : */
1509 : void constrain_element_residual (DenseVector<Number> & rhs,
1510 : std::vector<dof_id_type> & elem_dofs,
1511 : NumericVector<Number> & solution_local) const;
1512 :
1513 : /**
1514 : * Constrains a dyadic element matrix B = v w'. This method
1515 : * requires the element matrix to be square, in which case the
1516 : * elem_dofs correspond to the global DOF indices of both the rows
1517 : * and columns of the element matrix. For this case the rows and
1518 : * columns of the matrix necessarily correspond to variables of the
1519 : * same approximation order.
1520 : */
1521 : void constrain_element_dyad_matrix (DenseVector<Number> & v,
1522 : DenseVector<Number> & w,
1523 : std::vector<dof_id_type> & row_dofs,
1524 : bool asymmetric_constraint_rows = true) const;
1525 :
1526 : /**
1527 : * Does not actually constrain anything, but modifies \p dofs in the
1528 : * same way as any of the constrain functions would do, i.e. adds
1529 : * those dofs in terms of which any of the existing dofs is
1530 : * constrained.
1531 : */
1532 : void constrain_nothing (std::vector<dof_id_type> & dofs) const;
1533 :
1534 : /**
1535 : * Constrains the numeric vector \p v, which represents a solution defined on
1536 : * the mesh. This may need to be used after a linear solve, if your linear
1537 : * solver's solutions do not satisfy your DoF constraints to a tight enough
1538 : * tolerance.
1539 : *
1540 : * If \p v == nullptr, the system solution vector is constrained
1541 : *
1542 : * If \p homogeneous == true, heterogeneous constraints are enforced
1543 : * as if they were homogeneous. This might be appropriate for e.g. a
1544 : * vector representing a difference between two
1545 : * heterogeneously-constrained solutions.
1546 : */
1547 : void enforce_constraints_exactly (const System & system,
1548 : NumericVector<Number> * v = nullptr,
1549 : bool homogeneous = false) const;
1550 :
1551 : /**
1552 : * Heterogeneously constrains the numeric vector \p v, which
1553 : * represents an adjoint solution defined on the mesh for quantity
1554 : * fo interest \p q. For homogeneous constraints, use \p
1555 : * enforce_constraints_exactly instead
1556 : */
1557 : void enforce_adjoint_constraints_exactly (NumericVector<Number> & v,
1558 : unsigned int q) const;
1559 :
1560 : void enforce_constraints_on_residual (const NonlinearImplicitSystem & system,
1561 : NumericVector<Number> * rhs,
1562 : NumericVector<Number> const * solution,
1563 : bool homogeneous = true) const;
1564 :
1565 : void enforce_constraints_on_jacobian (const NonlinearImplicitSystem & system,
1566 : SparseMatrix<Number> * jac) const;
1567 :
1568 : #ifdef LIBMESH_ENABLE_PERIODIC
1569 :
1570 : //--------------------------------------------------------------------
1571 : // PeriodicBoundary-specific methods
1572 :
1573 : /**
1574 : * Adds a copy of the specified periodic boundary to the system.
1575 : */
1576 : void add_periodic_boundary (const PeriodicBoundaryBase & periodic_boundary);
1577 :
1578 : /**
1579 : * Add a periodic boundary pair
1580 : *
1581 : * \param boundary - primary boundary
1582 : * \param inverse_boundary - inverse boundary
1583 : */
1584 : void add_periodic_boundary (const PeriodicBoundaryBase & boundary, const PeriodicBoundaryBase & inverse_boundary);
1585 :
1586 : /**
1587 : * \returns \p true if the boundary given by \p boundaryid is periodic,
1588 : * false otherwise
1589 : */
1590 : bool is_periodic_boundary (const boundary_id_type boundaryid) const;
1591 :
1592 : PeriodicBoundaries * get_periodic_boundaries()
1593 : {
1594 : return _periodic_boundaries.get();
1595 : }
1596 :
1597 : const PeriodicBoundaries * get_periodic_boundaries() const
1598 : {
1599 : return _periodic_boundaries.get();
1600 : }
1601 :
1602 : #endif // LIBMESH_ENABLE_PERIODIC
1603 :
1604 :
1605 : #ifdef LIBMESH_ENABLE_DIRICHLET
1606 :
1607 : //--------------------------------------------------------------------
1608 : // DirichletBoundary-specific methods
1609 :
1610 : /**
1611 : * Adds a copy of the specified Dirichlet boundary to the system.
1612 : *
1613 : * The constraints implied by DirichletBoundary objects are imposed
1614 : * in the same order in which DirichletBoundary objects are added to
1615 : * the DofMap. When multiple DirichletBoundary objects would impose
1616 : * competing constraints on a given DOF, the *first*
1617 : * DirichletBoundary to constrain the DOF "wins". This distinction
1618 : * is important when e.g. two surfaces (sidesets) intersect. The
1619 : * nodes on the intersection will be constrained according to
1620 : * whichever sideset's DirichletBoundary object was added to the
1621 : * DofMap first.
1622 : */
1623 : void add_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary);
1624 :
1625 : /**
1626 : * Adds a copy of the specified Dirichlet boundary to the system,
1627 : * corresponding to the adjoint problem defined by Quantity of
1628 : * Interest \p q.
1629 : */
1630 : void add_adjoint_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary,
1631 : unsigned int q);
1632 :
1633 : /**
1634 : * Removes the specified Dirichlet boundary from the system.
1635 : */
1636 : void remove_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary);
1637 :
1638 : /**
1639 : * Removes from the system the specified Dirichlet boundary for the
1640 : * adjoint equation defined by Quantity of interest index q
1641 : */
1642 : void remove_adjoint_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary,
1643 : unsigned int q);
1644 :
1645 : const DirichletBoundaries * get_dirichlet_boundaries() const
1646 : {
1647 : return _dirichlet_boundaries.get();
1648 : }
1649 :
1650 24 : DirichletBoundaries * get_dirichlet_boundaries()
1651 : {
1652 24 : return _dirichlet_boundaries.get();
1653 : }
1654 :
1655 : bool has_adjoint_dirichlet_boundaries(unsigned int q) const;
1656 :
1657 : const DirichletBoundaries *
1658 : get_adjoint_dirichlet_boundaries(unsigned int q) const;
1659 :
1660 : DirichletBoundaries *
1661 : get_adjoint_dirichlet_boundaries(unsigned int q);
1662 :
1663 : /**
1664 : * Check that all the ids in dirichlet_bcids are actually present in the mesh.
1665 : * If not, this will throw an error.
1666 : */
1667 : void check_dirichlet_bcid_consistency (const MeshBase & mesh,
1668 : const DirichletBoundary & boundary) const;
1669 : #endif // LIBMESH_ENABLE_DIRICHLET
1670 :
1671 :
1672 : #ifdef LIBMESH_ENABLE_AMR
1673 :
1674 : //--------------------------------------------------------------------
1675 : // AMR-specific methods
1676 :
1677 : /**
1678 : * After a mesh is refined and repartitioned it is possible that the
1679 : * \p _send_list will need to be augmented. This is the case when an
1680 : * element is refined and its children end up on different processors
1681 : * than the parent. These children will need values from the parent
1682 : * when projecting the solution onto the refined mesh, hence the parent's
1683 : * DOF indices need to be included in the \p _send_list.
1684 : */
1685 : // void augment_send_list_for_projection(const MeshBase &);
1686 :
1687 : #ifdef LIBMESH_ENABLE_AMR
1688 :
1689 : /**
1690 : * Fills the vector di with the global degree of freedom indices
1691 : * for the element using the \p DofMap::old_dof_object.
1692 : * If no variable number is specified then all
1693 : * variables are returned.
1694 : */
1695 : void old_dof_indices (const Elem * const elem,
1696 : std::vector<dof_id_type> & di,
1697 : const unsigned int vn = libMesh::invalid_uint) const;
1698 :
1699 : #endif // LIBMESH_ENABLE_AMR
1700 :
1701 : /**
1702 : * Constrains degrees of freedom on side \p s of element \p elem which
1703 : * correspond to variable number \p var and to p refinement levels
1704 : * above \p p.
1705 : */
1706 : void constrain_p_dofs (unsigned int var,
1707 : const Elem * elem,
1708 : unsigned int s,
1709 : unsigned int p);
1710 :
1711 : #endif // LIBMESH_ENABLE_AMR
1712 :
1713 : /**
1714 : * Reinitialize the underlying data structures conformal to the current mesh.
1715 : */
1716 : void reinit
1717 : (MeshBase & mesh,
1718 : const std::map<const Node *, std::set<subdomain_id_type>> &
1719 : constraining_subdomains);
1720 :
1721 : /**
1722 : * Free all new memory associated with the object, but restore its
1723 : * original state, with the mesh pointer and any default ghosting.
1724 : */
1725 : virtual void clear () override;
1726 :
1727 : /**
1728 : * Prints summary info about the sparsity bandwidth and constraints.
1729 : */
1730 : void print_info(std::ostream & os=libMesh::out) const;
1731 :
1732 : /**
1733 : * Gets summary info about the sparsity bandwidth and constraints.
1734 : */
1735 : std::string get_info() const;
1736 :
1737 : /**
1738 : * Degree of freedom coupling. If left empty each DOF
1739 : * couples to all others. Can be used to reduce memory
1740 : * requirements for sparse matrices. DOF 0 might only
1741 : * couple to itself, in which case \p dof_coupling(0,0)
1742 : * should be 1 and \p dof_coupling(0,j) = 0 for j not equal
1743 : * to 0.
1744 : *
1745 : * This variable is named as though it were class private,
1746 : * but it is in the public interface. Also there are no
1747 : * public methods for accessing it... This typically means
1748 : * you should only use it if you know what you are doing.
1749 : */
1750 : CouplingMatrix * _dof_coupling;
1751 :
1752 : /**
1753 : * \returns The number of the system we are responsible for.
1754 : */
1755 : unsigned int sys_number() const;
1756 :
1757 : /**
1758 : * Builds a sparsity pattern for matrices using the current
1759 : * degree-of-freedom numbering and coupling.
1760 : *
1761 : * By default, ignores constraint equations, for build speed; this
1762 : * is valid for the combination of !need_full_sparsity_pattern and
1763 : * constraints which only come from periodic boundary conditions and
1764 : * adaptive mesh refinement, where matrix constraint adds some
1765 : * matrix entries but removes equally many (or more) other entries.
1766 : *
1767 : * Can be told to calculate sparsity for the constrained matrix,
1768 : * which may be necessary in the case of spline control node
1769 : * constraints or sufficiently many user constraints.
1770 : *
1771 : * If \p reduced is provided, the pattern couples only the degrees
1772 : * of freedom that \p reduced retains, giving the nonzero counts of
1773 : * matrices in that reduced space.
1774 : */
1775 : std::unique_ptr<SparsityPattern::Build> build_sparsity(const MeshBase & mesh,
1776 : bool calculate_constrained = false,
1777 : const ReducedDofMap * reduced = nullptr) const;
1778 :
1779 : /**
1780 : * Set whether the given variable group should be p-refined on a
1781 : * p-refined Elem. This changes the FEType of the variable group to
1782 : * enable or disable p-refinement.
1783 : */
1784 : void should_p_refine(unsigned int g, bool p_refine);
1785 :
1786 : /**
1787 : * Whether the given variable group should be p-refined
1788 : */
1789 : bool should_p_refine(unsigned int g) const;
1790 :
1791 : /**
1792 : * Whether the given variable should be p-refined
1793 : */
1794 : bool should_p_refine_var(unsigned int var) const;
1795 :
1796 : // Prevent bad user implicit conversions
1797 : void should_p_refine(FEFamily, bool) = delete;
1798 : void should_p_refine(Order, bool) = delete;
1799 : bool should_p_refine(FEFamily) const = delete;
1800 : bool should_p_refine(Order) const = delete;
1801 :
1802 : /**
1803 : * Add a static condensation class
1804 : */
1805 : void create_static_condensation(MeshBase & mesh, System & system);
1806 :
1807 : /**
1808 : * Checks whether we have static condensation
1809 : */
1810 9716 : bool has_static_condensation() const { return _sc.get(); }
1811 :
1812 : /**
1813 : * @returns the static condensation class. This should have been already added with a call to \p
1814 : * add_static_condensation()
1815 : */
1816 : StaticCondensationDofMap & get_static_condensation();
1817 :
1818 : /**
1819 : * @returns the static condensation class. This should have been already added with a call to \p
1820 : * add_static_condensation()
1821 : */
1822 : const StaticCondensationDofMap & get_static_condensation() const;
1823 :
1824 : /**
1825 : * Calls reinit on the static condensation map if it exists
1826 : */
1827 : void reinit_static_condensation();
1828 :
1829 : /**
1830 : * Add a map of the unconstrained degrees of freedom, so that solves can eliminate the
1831 : * constrained ones
1832 : */
1833 : void create_constraint_elimination(MeshBase & mesh, System & system);
1834 :
1835 : /**
1836 : * Checks whether we eliminate constrained degrees of freedom
1837 : */
1838 3132 : bool has_constraint_elimination() const { return _constraint_elimination.get(); }
1839 :
1840 : /**
1841 : * @returns The map of the unconstrained degrees of freedom. This should have been already added
1842 : * with a call to \p create_constraint_elimination()
1843 : */
1844 : const ConstraintEliminationDofMap & get_constraint_elimination() const;
1845 :
1846 : /**
1847 : * Calls reinit on the constraint elimination map if it exists
1848 : */
1849 : void reinit_constraint_elimination();
1850 :
1851 : private:
1852 :
1853 : /**
1854 : * Retrieve the array variable bounds for a given variable \p vi. This variable may
1855 : * lie anywhere within an array variable range. An 'array variable' is simply a sequence
1856 : * of contiguous variable numbers defined by pair where the first member of the pair
1857 : * is the first number in the variable sequence and the second member of the pair is
1858 : * the number of the last variable in the sequence plus one. Array variables may be
1859 : * used in tandem with variable grouping by downstream code to build optimized physics
1860 : * kernels since each variable in the array will have the same shape functions.
1861 : *
1862 : * We note that we store array variables as a container of the above described pairs. Within
1863 : * this API we will do a binary search such that the complexity is O(log(N)) where N is the
1864 : * number of array variables present in \p this
1865 : */
1866 : const std::pair<unsigned int, unsigned int> &
1867 : get_variable_array(unsigned int vi) const;
1868 :
1869 : /**
1870 : * Helper function that gets the dof indices on the current element
1871 : * for a non-SCALAR type variable, where the variable is identified
1872 : * by its variable group number \p vg and its offset \p vig from the
1873 : * first variable in that group.
1874 : *
1875 : * In DEBUG mode, the tot_size parameter will add up the total
1876 : * number of dof indices that should have been added to di, and v
1877 : * will be the variable number corresponding to vg and vig.
1878 : */
1879 : void _dof_indices (const Elem & elem,
1880 : int p_level,
1881 : std::vector<dof_id_type> & di,
1882 : const unsigned int vg,
1883 : const unsigned int vig,
1884 : const Node * const * nodes,
1885 : unsigned int n_nodes,
1886 : const unsigned int v
1887 : #ifdef DEBUG
1888 : ,
1889 : std::size_t & tot_size
1890 : #endif
1891 : ) const;
1892 :
1893 : /**
1894 : * As above except a \p field_dofs_functor must be provided. This method is useful when the caller
1895 : * wants to do more than simply fill a degree of freedom container
1896 : * @param field_dofs_functor This functor has the interface:
1897 : * void field_dofs_functor(const Elem & elem,
1898 : * const unsigned int node_num,
1899 : * const unsigned int var_num,
1900 : * std::vector<dof_id_type> & di,
1901 : * const dof_id_type field_dof)
1902 : * where \p field_dof represents a field degree of freedom to act on and
1903 : * is associated with \p node_num and \p var_num. If the degree of
1904 : * freedom is elemental than \p node_num will be \p invalid_uint. \p di
1905 : * is the degree of freedom container provided to the \p _dof_indices
1906 : * method
1907 : */
1908 : template <typename FieldDofsFunctor>
1909 : void _dof_indices (const Elem & elem,
1910 : int p_level,
1911 : std::vector<dof_id_type> & di,
1912 : const unsigned int vg,
1913 : const unsigned int vig,
1914 : const Node * const * nodes,
1915 : unsigned int n_nodes,
1916 : const unsigned int v,
1917 : #ifdef DEBUG
1918 : std::size_t & tot_size,
1919 : #endif
1920 : FieldDofsFunctor field_dofs_functor) const;
1921 :
1922 : /**
1923 : * Helper function that implements the element-nodal versions of
1924 : * dof_indices and old_dof_indices
1925 : */
1926 : void _node_dof_indices (const Elem & elem,
1927 : unsigned int n,
1928 : const DofObject & obj,
1929 : std::vector<dof_id_type> & di,
1930 : const unsigned int vn) const;
1931 :
1932 : /**
1933 : * Invalidates all active DofObject dofs for this system
1934 : */
1935 : void invalidate_dofs(MeshBase & mesh) const;
1936 :
1937 : /**
1938 : * \returns The Node pointer with index \p i from the \p mesh.
1939 : */
1940 : DofObject * node_ptr(MeshBase & mesh, dof_id_type i) const;
1941 :
1942 : /**
1943 : * \returns The Elem pointer with index \p i from the \p mesh.
1944 : */
1945 : DofObject * elem_ptr(MeshBase & mesh, dof_id_type i) const;
1946 :
1947 : /**
1948 : * A member function type like \p node_ptr() or \p elem_ptr().
1949 : */
1950 : typedef DofObject * (DofMap::*dofobject_accessor)
1951 : (MeshBase & mesh, dof_id_type i) const;
1952 :
1953 : /**
1954 : * Helper function for distributing dofs in parallel
1955 : */
1956 : template<typename iterator_type>
1957 : void set_nonlocal_dof_objects(iterator_type objects_begin,
1958 : iterator_type objects_end,
1959 : MeshBase & mesh,
1960 : dofobject_accessor objects);
1961 :
1962 : /**
1963 : * We may have mesh constraint rows with dependent nodes in one
1964 : * subdomain but dependency nodes in another subdomain, and we may
1965 : * have variables whose subdomain restriction includes the dependent
1966 : * subdomain but not the dependency. In those cases we need to
1967 : * place degrees of freedom on dependency nodes anyway.
1968 : *
1969 : * The set value for node n will include all subdomain ids of
1970 : * elements with nodes in subdomains constrained by n.
1971 : *
1972 : * We use a map<set> rather than a multimap here because we expect
1973 : * to be inserting the same subdomain multiple times and we don't
1974 : * need duplicate values.
1975 : */
1976 : std::map<const Node *, std::set<subdomain_id_type>>
1977 : calculate_constraining_subdomains();
1978 :
1979 : /**
1980 : * Distributes the global degrees of freedom, for dofs on
1981 : * this processor. In this format the local
1982 : * degrees of freedom are in a contiguous block for each
1983 : * variable in the system.
1984 : * Starts at index next_free_dof, and increments it to
1985 : * the post-final index.
1986 : *
1987 : * Uses the provided constraining_subdomains map from
1988 : * calculate_constraining_subdomains() to ensure allocation of all
1989 : * DoFs on constraining nodes.
1990 : */
1991 : void distribute_local_dofs_var_major
1992 : (dof_id_type & next_free_dof,
1993 : MeshBase & mesh,
1994 : const std::map<const Node *, std::set<subdomain_id_type>> &
1995 : constraining_subdomains);
1996 :
1997 : /**
1998 : * Distributes the global degrees of freedom for dofs on this
1999 : * processor. In this format all the degrees of freedom at a
2000 : * node/element are in contiguous blocks. Starts at index \p
2001 : * next_free_dof, and increments it to the post-final index. If \p
2002 : * build_send_list is \p true, builds the send list. If \p false,
2003 : * clears and reserves the send list.
2004 : *
2005 : * Uses the provided constraining_subdomains map from
2006 : * calculate_constraining_subdomains() to ensure allocation of all
2007 : * DoFs on constraining nodes.
2008 : *
2009 : * \note The degrees of freedom for a given variable are not in
2010 : * contiguous blocks, as in the case of \p distribute_local_dofs_var_major.
2011 : */
2012 : void distribute_local_dofs_node_major
2013 : (dof_id_type & next_free_dof,
2014 : MeshBase & mesh,
2015 : const std::map<const Node *, std::set<subdomain_id_type>> &
2016 : constraining_subdomains);
2017 :
2018 : /*
2019 : * Helper method for the above two to count + distriubte SCALAR dofs
2020 : */
2021 : void distribute_scalar_dofs (dof_id_type & next_free_dof);
2022 :
2023 : #ifdef DEBUG
2024 : /*
2025 : * Internal assertions for distribute_local_dofs_*
2026 : */
2027 : void assert_no_nodes_missed(MeshBase & mesh);
2028 : #endif
2029 :
2030 : /*
2031 : * A utility method for obtaining a set of elements to ghost along
2032 : * with merged coupling matrices.
2033 : */
2034 : typedef std::set<std::unique_ptr<CouplingMatrix>, Utility::CompareUnderlying> CouplingMatricesSet;
2035 : static void
2036 : merge_ghost_functor_outputs (GhostingFunctor::map_type & elements_to_ghost,
2037 : CouplingMatricesSet & temporary_coupling_matrices,
2038 : const GhostingFunctorIterator & gf_begin,
2039 : const GhostingFunctorIterator & gf_end,
2040 : const MeshBase::const_element_iterator & elems_begin,
2041 : const MeshBase::const_element_iterator & elems_end,
2042 : processor_id_type p);
2043 :
2044 : /**
2045 : * Adds entries to the \p _send_list vector corresponding to DoFs
2046 : * on elements neighboring the current processor.
2047 : */
2048 : void add_neighbors_to_send_list(MeshBase & mesh);
2049 :
2050 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
2051 :
2052 : /**
2053 : * Build the constraint matrix C associated with the element
2054 : * degree of freedom indices elem_dofs. The optional parameter
2055 : * \p called_recursively should be left at the default value
2056 : * \p false. This is used to handle the special case of
2057 : * an element's degrees of freedom being constrained in terms
2058 : * of other, local degrees of freedom. The usual case is
2059 : * for an elements DOFs to be constrained by some other,
2060 : * external DOFs.
2061 : */
2062 : void build_constraint_matrix (DenseMatrix<Number> & C,
2063 : std::vector<dof_id_type> & elem_dofs,
2064 : const bool called_recursively=false) const;
2065 :
2066 : /**
2067 : * Build the constraint matrix C and the forcing vector H
2068 : * associated with the element degree of freedom indices elem_dofs.
2069 : * The optional parameter \p called_recursively should be left at
2070 : * the default value \p false. This is used to handle the special
2071 : * case of an element's degrees of freedom being constrained in
2072 : * terms of other, local degrees of freedom. The usual case is for
2073 : * an elements DOFs to be constrained by some other, external DOFs
2074 : * and/or Dirichlet conditions.
2075 : *
2076 : * The forcing vector will depend on which solution's heterogeneous
2077 : * constraints are being applied. For the default \p qoi_index this
2078 : * will be the primal solution; for \p qoi_index >= 0 the
2079 : * corresponding adjoint solution's constraints will be used.
2080 : */
2081 : void build_constraint_matrix_and_vector (DenseMatrix<Number> & C,
2082 : DenseVector<Number> & H,
2083 : std::vector<dof_id_type> & elem_dofs,
2084 : int qoi_index = -1,
2085 : const bool called_recursively=false) const;
2086 :
2087 : /**
2088 : * Finds all the DOFS associated with the element DOFs elem_dofs.
2089 : * This will account for off-element couplings via hanging nodes.
2090 : */
2091 : void find_connected_dofs (std::vector<dof_id_type> & elem_dofs) const;
2092 :
2093 : /**
2094 : * Finds all the DofObjects associated with the set in \p objs.
2095 : * This will account for off-element couplings via hanging nodes.
2096 : */
2097 : void find_connected_dof_objects (std::vector<const DofObject *> & objs) const;
2098 :
2099 : /**
2100 : * Adds entries to the \p _send_list vector corresponding to DoFs
2101 : * which are dependencies for constraint equations on the current
2102 : * processor.
2103 : */
2104 : void add_constraints_to_send_list(const MeshBase & mesh);
2105 :
2106 : /**
2107 : * Adds any spline constraints from the Mesh to our DoF constraints.
2108 : * If any Dirichlet constraints exist on spline-constrained nodes,
2109 : * l2-projects those constraints onto the spline basis.
2110 : */
2111 : void process_mesh_constraint_rows(const MeshBase & mesh);
2112 :
2113 : #endif // LIBMESH_ENABLE_CONSTRAINTS
2114 :
2115 : /**
2116 : * This flag indicates whether or not we do an opt-mode check for
2117 : * the presence of constraint loops, i.e. cases where the constraint
2118 : * graph is cyclic.
2119 : */
2120 : bool _error_on_constraint_loop;
2121 :
2122 : /**
2123 : * This flag indicates whether or not we explicitly take constraint
2124 : * equations into account when computing a sparsity pattern.
2125 : */
2126 : bool _constrained_sparsity_construction;
2127 :
2128 : /**
2129 : * The variables in this system/degree of freedom map
2130 : */
2131 : std::vector<Variable> _variables;
2132 :
2133 : /**
2134 : * The variable groups in this system/degree of freedom map
2135 : */
2136 : std::vector<VariableGroup> _variable_groups;
2137 :
2138 : /**
2139 : * The variable group number for each variable.
2140 : */
2141 : std::vector<unsigned int> _variable_group_numbers;
2142 :
2143 : /**
2144 : * A map from variable number to variable group number
2145 : */
2146 : std::unordered_map<unsigned int, unsigned int> _var_to_vg;
2147 :
2148 : /**
2149 : * The variable numbers corresponding to user-specified
2150 : * names, useful for name-based lookups.
2151 : */
2152 : std::map<std::string, unsigned int, std::less<>> _variable_numbers;
2153 :
2154 : /**
2155 : * Array variable information storage. For a given array "variable", the first member of the pair
2156 : * denotes the first variable number present in the array variable and the second member of the
2157 : * pair denotes the last variable number present in the array variables plus one
2158 : */
2159 : std::vector<std::pair<unsigned int, unsigned int>> _array_variables;
2160 :
2161 : /**
2162 : * \p true when \p VariableGroup structures should be automatically
2163 : * identified, \p false otherwise. Defaults to \p true.
2164 : */
2165 : bool _identify_variable_groups = true;
2166 :
2167 : /**
2168 : * The number of the system we manage DOFs for.
2169 : */
2170 : const unsigned int _sys_number;
2171 :
2172 : /**
2173 : * The mesh that system uses.
2174 : */
2175 : MeshBase & _mesh;
2176 :
2177 : /**
2178 : * Additional matrices handled by this object. These pointers do \e
2179 : * not handle the memory, instead, \p System, who
2180 : * told \p DofMap about them, owns them.
2181 : */
2182 : std::vector<SparseMatrix<Number> * > _matrices;
2183 :
2184 : /**
2185 : * First DOF index for SCALAR variable v, or garbage for non-SCALAR
2186 : * variable v
2187 : */
2188 : std::vector<dof_id_type> _first_scalar_df;
2189 :
2190 : /**
2191 : * A list containing all the global DOF indices that affect the
2192 : * solution on my processor.
2193 : */
2194 : std::vector<dof_id_type> _send_list;
2195 :
2196 : /**
2197 : * Function object to call to add extra entries to the sparsity pattern
2198 : */
2199 : SparsityPattern::AugmentSparsityPattern * _augment_sparsity_pattern;
2200 :
2201 : /**
2202 : * A function pointer to a function to call to add extra entries to the sparsity pattern
2203 : */
2204 : void (*_extra_sparsity_function)(SparsityPattern::Graph &,
2205 : std::vector<dof_id_type> & n_nz,
2206 : std::vector<dof_id_type> & n_oz,
2207 : void *);
2208 : /**
2209 : * A pointer associated with the extra sparsity that can optionally be passed in
2210 : */
2211 : void * _extra_sparsity_context;
2212 :
2213 : /**
2214 : * Function object to call to add extra entries to the send list
2215 : */
2216 : AugmentSendList * _augment_send_list;
2217 :
2218 : /**
2219 : * A function pointer to a function to call to add extra entries to the send list
2220 : */
2221 : void (*_extra_send_list_function)(std::vector<dof_id_type> &, void *);
2222 :
2223 : /**
2224 : * A pointer associated with the extra send list that can optionally be passed in
2225 : */
2226 : void * _extra_send_list_context;
2227 :
2228 : /**
2229 : * The default coupling GhostingFunctor, used to implement standard
2230 : * libMesh sparsity pattern construction.
2231 : *
2232 : * We use a std::unique_ptr here to reduce header dependencies.
2233 : */
2234 : std::unique_ptr<DefaultCoupling> _default_coupling;
2235 :
2236 : /**
2237 : * The default algebraic GhostingFunctor, used to implement standard
2238 : * libMesh send_list construction.
2239 : *
2240 : * We use a std::unique_ptr here to reduce header dependencies.
2241 : */
2242 : std::unique_ptr<DefaultCoupling> _default_evaluating;
2243 :
2244 : /**
2245 : * The list of all GhostingFunctor objects to be used when
2246 : * distributing ghosted vectors.
2247 : *
2248 : * The library should automatically refer these functors to the
2249 : * MeshBase, too, so any algebraically ghosted dofs will live on
2250 : * geometrically ghosted elements.
2251 : *
2252 : * Keep these in a vector so any parallel computation is done in the
2253 : * same order on all processors.
2254 : */
2255 : std::vector<GhostingFunctor *> _algebraic_ghosting_functors;
2256 :
2257 : /**
2258 : * The list of all GhostingFunctor objects to be used when
2259 : * coupling degrees of freedom in matrix sparsity patterns.
2260 : *
2261 : * These objects will *also* be used as algebraic ghosting functors,
2262 : * but not vice-versa.
2263 : *
2264 : * The library should automatically refer these functors to the
2265 : * MeshBase, too, so any dofs coupled to local dofs will live on
2266 : * geometrically ghosted elements.
2267 : */
2268 : std::vector<GhostingFunctor *> _coupling_functors;
2269 :
2270 : /**
2271 : * Hang on to references to any GhostingFunctor objects we were
2272 : * passed in shared_ptr form
2273 : */
2274 : std::map<GhostingFunctor *, std::shared_ptr<GhostingFunctor> > _shared_functors;
2275 :
2276 : /**
2277 : * Default false; set to true if any attached matrix requires a full
2278 : * sparsity pattern.
2279 : */
2280 : bool _need_full_sparsity_pattern;
2281 :
2282 : /**
2283 : * Default false; set to true if the dependencies of constrained ghost
2284 : * DOFs supported by local elements should also be ghosted
2285 : */
2286 : bool _need_ghost_constraints;
2287 :
2288 : /**
2289 : * The sparsity pattern of the global matrix. If
2290 : * need_full_sparsity_pattern is true, we save the entire sparse
2291 : * graph here. Otherwise we save just the n_nz and n_oz vectors.
2292 : */
2293 : std::unique_ptr<SparsityPattern::Build> _sp;
2294 :
2295 : /**
2296 : * The total number of SCALAR dofs associated to
2297 : * all SCALAR variables.
2298 : */
2299 : dof_id_type _n_SCALAR_dofs;
2300 :
2301 : #ifdef LIBMESH_ENABLE_AMR
2302 :
2303 : /**
2304 : * First old DOF index for SCALAR variable v, or garbage for
2305 : * non-SCALAR variable v
2306 : */
2307 : std::vector<dof_id_type> _first_old_scalar_df;
2308 : #endif
2309 :
2310 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
2311 : /**
2312 : * Data structure containing DOF constraints. The ith
2313 : * entry is the constraint matrix row for DOF i.
2314 : */
2315 : DofConstraints _dof_constraints, _stashed_dof_constraints;
2316 :
2317 : DofConstraintValueMap _primal_constraint_values;
2318 :
2319 : AdjointDofConstraintValues _adjoint_constraint_values;
2320 : #endif
2321 :
2322 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2323 : /**
2324 : * Data structure containing DofObject constraints.
2325 : */
2326 : NodeConstraints _node_constraints;
2327 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2328 :
2329 :
2330 : #ifdef LIBMESH_ENABLE_PERIODIC
2331 : /**
2332 : * Data structure containing periodic boundaries. The ith
2333 : * entry is the constraint matrix row for boundaryid i.
2334 : */
2335 : std::unique_ptr<PeriodicBoundaries> _periodic_boundaries;
2336 : #endif
2337 :
2338 : #ifdef LIBMESH_ENABLE_DIRICHLET
2339 : /**
2340 : * Data structure containing Dirichlet functions. The ith
2341 : * entry is the constraint matrix row for boundaryid i.
2342 : */
2343 : std::unique_ptr<DirichletBoundaries> _dirichlet_boundaries;
2344 :
2345 : /**
2346 : * Data structure containing Dirichlet functions. The ith
2347 : * entry is the constraint matrix row for boundaryid i.
2348 : */
2349 : std::vector<std::unique_ptr<DirichletBoundaries>> _adjoint_dirichlet_boundaries;
2350 : #endif
2351 :
2352 : friend class SparsityPattern::Build;
2353 :
2354 : /**
2355 : * Bools to indicate if we override the --implicit_neighbor_dofs
2356 : * commandline options.
2357 : */
2358 : bool _implicit_neighbor_dofs_initialized;
2359 : bool _implicit_neighbor_dofs;
2360 :
2361 : /**
2362 : * Flag which determines whether we should do some additional
2363 : * checking of the consistency of the DirichletBoundary objects
2364 : * added by the user. Defaults to true, but can be disabled in cases
2365 : * where you only want to add DirichletBoundary objects "locally"
2366 : * and can guarantee that no repartitioning will be done, since
2367 : * repartitioning could cause processors to own new boundary sides
2368 : * for which they no longer have the proper DirichletBoundary
2369 : * objects stored.
2370 : */
2371 : bool _verify_dirichlet_bc_consistency;
2372 :
2373 : /// Static condensation class
2374 : std::unique_ptr<StaticCondensationDofMap> _sc;
2375 :
2376 : /// The map of the unconstrained degrees of freedom, when solves eliminate constrained ones
2377 : std::unique_ptr<ConstraintEliminationDofMap> _constraint_elimination;
2378 : };
2379 :
2380 :
2381 : // ------------------------------------------------------------
2382 : // Dof Map inline member functions
2383 : inline
2384 44489061 : unsigned int DofMap::sys_number() const
2385 : {
2386 73973392 : return _sys_number;
2387 : }
2388 :
2389 :
2390 :
2391 : inline
2392 48981621 : const VariableGroup & DofMap::variable_group (const unsigned int g) const
2393 : {
2394 48981621 : libmesh_assert_less (g, _variable_groups.size());
2395 :
2396 98132974 : return _variable_groups[g];
2397 : }
2398 :
2399 :
2400 :
2401 : inline
2402 21452198 : const Variable & DofMap::variable (const unsigned int c) const
2403 : {
2404 6531646 : libmesh_assert_less (c, _variables.size());
2405 :
2406 26768587 : return _variables[c];
2407 : }
2408 :
2409 :
2410 :
2411 : inline
2412 : Order DofMap::variable_order (const unsigned int c) const
2413 : {
2414 : libmesh_assert_less (c, _variables.size());
2415 :
2416 : return _variables[c].type().order;
2417 : }
2418 :
2419 :
2420 :
2421 : inline
2422 : Order DofMap::variable_group_order (const unsigned int vg) const
2423 : {
2424 : libmesh_assert_less (vg, _variable_groups.size());
2425 :
2426 : return _variable_groups[vg].type().order;
2427 : }
2428 :
2429 :
2430 :
2431 : inline
2432 1883052 : const FEType & DofMap::variable_type (const unsigned int c) const
2433 : {
2434 1883052 : libmesh_assert_less (c, _variables.size());
2435 :
2436 5855027 : return _variables[c].type();
2437 : }
2438 :
2439 :
2440 :
2441 : inline
2442 : const FEType & DofMap::variable_group_type (const unsigned int vg) const
2443 : {
2444 : libmesh_assert_less (vg, _variable_groups.size());
2445 :
2446 : return _variable_groups[vg].type();
2447 : }
2448 :
2449 :
2450 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
2451 :
2452 :
2453 : inline
2454 1427761 : bool DofMap::is_constrained_node (const Node *
2455 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2456 : node
2457 : #endif
2458 : ) const
2459 : {
2460 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2461 1427761 : if (_node_constraints.count(node))
2462 17700 : return true;
2463 : #endif
2464 :
2465 1410061 : return false;
2466 : }
2467 :
2468 :
2469 : inline
2470 52849628 : bool DofMap::is_constrained_dof (const dof_id_type dof) const
2471 : {
2472 52849628 : if (_dof_constraints.count(dof))
2473 5774993 : return true;
2474 :
2475 47075320 : return false;
2476 : }
2477 :
2478 :
2479 : inline
2480 6186 : bool DofMap::has_heterogeneous_adjoint_constraints (const unsigned int qoi_num) const
2481 : {
2482 : AdjointDofConstraintValues::const_iterator it =
2483 3072 : _adjoint_constraint_values.find(qoi_num);
2484 9258 : if (it == _adjoint_constraint_values.end())
2485 1618 : return false;
2486 4368 : if (it->second.empty())
2487 2880 : return false;
2488 :
2489 14 : return true;
2490 : }
2491 :
2492 :
2493 : inline
2494 1244328 : Number DofMap::has_heterogeneous_adjoint_constraint (const unsigned int qoi_num,
2495 : const dof_id_type dof) const
2496 : {
2497 : AdjointDofConstraintValues::const_iterator it =
2498 1093446 : _adjoint_constraint_values.find(qoi_num);
2499 1244328 : if (it != _adjoint_constraint_values.end())
2500 : {
2501 : DofConstraintValueMap::const_iterator rhsit =
2502 104428 : it->second.find(dof);
2503 255310 : if (rhsit == it->second.end())
2504 154302 : return 0;
2505 : else
2506 700 : return rhsit->second;
2507 : }
2508 :
2509 989018 : return 0;
2510 : }
2511 :
2512 :
2513 :
2514 : inline
2515 : DofConstraintValueMap & DofMap::get_primal_constraint_values()
2516 : {
2517 : return _primal_constraint_values;
2518 : }
2519 :
2520 :
2521 :
2522 : #else
2523 :
2524 : //--------------------------------------------------------------------
2525 : // Constraint-specific methods get inlined into nothing if
2526 : // constraints are disabled, so there's no reason for users not to
2527 : // use them.
2528 :
2529 : inline void DofMap::constrain_element_matrix (DenseMatrix<Number> &,
2530 : std::vector<dof_id_type> &,
2531 : bool) const {}
2532 :
2533 : inline void DofMap::constrain_element_matrix (DenseMatrix<Number> &,
2534 : std::vector<dof_id_type> &,
2535 : std::vector<dof_id_type> &,
2536 : bool) const {}
2537 :
2538 : inline void DofMap::constrain_element_vector (DenseVector<Number> &,
2539 : std::vector<dof_id_type> &,
2540 : bool) const {}
2541 :
2542 : inline void DofMap::constrain_element_matrix_and_vector (DenseMatrix<Number> &,
2543 : DenseVector<Number> &,
2544 : std::vector<dof_id_type> &,
2545 : bool) const {}
2546 :
2547 : inline void DofMap::heterogeneously_constrain_element_matrix_and_vector
2548 : (DenseMatrix<Number> &, DenseVector<Number> &,
2549 : std::vector<dof_id_type> &, bool, int) const {}
2550 :
2551 : inline void DofMap::heterogeneously_constrain_element_vector
2552 : (const DenseMatrix<Number> &, DenseVector<Number> &,
2553 : std::vector<dof_id_type> &, bool, int) const {}
2554 :
2555 : inline void DofMap::constrain_element_dyad_matrix (DenseVector<Number> &,
2556 : DenseVector<Number> &,
2557 : std::vector<dof_id_type> &,
2558 : bool) const {}
2559 :
2560 : inline void DofMap::constrain_nothing (std::vector<dof_id_type> &) const {}
2561 :
2562 : inline void DofMap::enforce_constraints_exactly (const System &,
2563 : NumericVector<Number> *,
2564 : bool) const {}
2565 :
2566 : inline void DofMap::enforce_adjoint_constraints_exactly (NumericVector<Number> &,
2567 : unsigned int) const {}
2568 :
2569 :
2570 : inline void DofMap::enforce_constraints_on_residual
2571 : (const NonlinearImplicitSystem &,
2572 : NumericVector<Number> *,
2573 : NumericVector<Number> const *,
2574 : bool) const {}
2575 :
2576 : inline void DofMap::enforce_constraints_on_jacobian
2577 : (const NonlinearImplicitSystem &,
2578 : SparseMatrix<Number> *) const {}
2579 :
2580 : #endif // LIBMESH_ENABLE_CONSTRAINTS
2581 :
2582 :
2583 :
2584 : inline
2585 : void DofMap::set_constrained_sparsity_construction(bool use_constraints)
2586 : {
2587 : // This got only partly finished...
2588 : if (use_constraints)
2589 : libmesh_not_implemented();
2590 :
2591 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
2592 : _constrained_sparsity_construction = use_constraints;
2593 : #endif
2594 : libmesh_ignore(use_constraints);
2595 : }
2596 :
2597 : inline
2598 : void DofMap::full_sparsity_pattern_needed()
2599 : {
2600 : _need_full_sparsity_pattern = true;
2601 : }
2602 :
2603 : inline
2604 : void DofMap::ghost_constraints_needed()
2605 : {
2606 : _need_ghost_constraints = true;
2607 : }
2608 :
2609 : inline
2610 : bool DofMap::constrained_sparsity_construction()
2611 : {
2612 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
2613 : return _constrained_sparsity_construction;
2614 : #else
2615 : return true;
2616 : #endif
2617 : }
2618 :
2619 : inline
2620 : void DofMap::should_p_refine(const unsigned int g, const bool p_refine)
2621 : {
2622 : #ifdef LIBMESH_ENABLE_AMR
2623 : VariableGroup & var = _variable_groups[g];
2624 : var.type().p_refinement = p_refine;
2625 :
2626 : for (auto v : make_range(var.first_scalar_number(0),
2627 : var.first_scalar_number(0) +
2628 : var.n_variables()))
2629 : this->_variables[v].type().p_refinement = p_refine;
2630 :
2631 :
2632 : #else
2633 : libmesh_ignore(g, p_refine);
2634 : #endif
2635 : }
2636 :
2637 : inline
2638 : bool DofMap::should_p_refine(const unsigned int g) const
2639 : {
2640 : #ifdef LIBMESH_ENABLE_AMR
2641 : const VariableGroup & var = this->variable_group(g);
2642 : return var.type().p_refinement;
2643 : #else
2644 : libmesh_ignore(g);
2645 : return false;
2646 : #endif
2647 : }
2648 :
2649 : inline
2650 : unsigned int DofMap::var_group_from_var_number(const unsigned int var_num) const
2651 : {
2652 : libmesh_assert(var_num < n_variables());
2653 : return libmesh_map_find(_var_to_vg, var_num);
2654 : }
2655 :
2656 : inline
2657 : bool DofMap::should_p_refine_var(const unsigned int var) const
2658 : {
2659 : #ifdef LIBMESH_ENABLE_AMR
2660 : const auto vg = this->var_group_from_var_number(var);
2661 : return this->should_p_refine(vg);
2662 : #else
2663 : libmesh_ignore(var);
2664 : return false;
2665 : #endif
2666 : }
2667 :
2668 : template <typename FieldDofsFunctor>
2669 58265858 : void DofMap::_dof_indices (const Elem & elem,
2670 : int p_level,
2671 : std::vector<dof_id_type> & di,
2672 : const unsigned int vg,
2673 : const unsigned int vig,
2674 : const Node * const * nodes,
2675 : unsigned int n_nodes,
2676 : const unsigned int v,
2677 : #ifdef DEBUG
2678 : std::size_t & tot_size,
2679 : #endif
2680 : FieldDofsFunctor field_dofs_functor) const
2681 : {
2682 24019846 : const VariableGroup & var = this->variable_group(vg);
2683 :
2684 58265858 : if (var.active_on_subdomain(elem.subdomain_id()))
2685 : {
2686 58238830 : const ElemType type = elem.type();
2687 24072795 : const unsigned int sys_num = this->sys_number();
2688 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2689 31610937 : const bool is_inf = elem.infinite();
2690 : #endif
2691 :
2692 : const bool extra_hanging_dofs =
2693 58238830 : FEInterface::extra_hanging_dofs(var.type());
2694 :
2695 58238830 : FEType fe_type = var.type();
2696 :
2697 58238830 : const bool add_p_level = fe_type.p_refinement;
2698 :
2699 : #ifdef DEBUG
2700 : // The number of dofs per element is non-static for subdivision FE
2701 24009473 : if (var.type().family == SUBDIVISION)
2702 3936 : tot_size += n_nodes;
2703 : else
2704 : // FIXME: Is the passed-in p_level just elem.p_level()? If so,
2705 : // this seems redundant.
2706 24005537 : tot_size += FEInterface::n_dofs(fe_type, add_p_level*p_level, &elem);
2707 : #endif
2708 :
2709 : // The total Order is not required when getting the function
2710 : // pointer, it is only needed when the function is called (see
2711 : // below).
2712 : const FEInterface::n_dofs_at_node_ptr ndan =
2713 58238830 : FEInterface::n_dofs_at_node_function(fe_type, &elem);
2714 :
2715 : // Get the node-based DOF numbers
2716 366846806 : for (unsigned int n=0; n != n_nodes; n++)
2717 : {
2718 308607976 : const Node & node = *nodes[n];
2719 :
2720 : // Cache the intermediate lookups that are common to every
2721 : // component
2722 : #ifdef DEBUG
2723 : const std::pair<unsigned int, unsigned int>
2724 121955875 : vg_and_offset = node.var_to_vg_and_offset(sys_num,v);
2725 121955875 : libmesh_assert_equal_to (vg, vg_and_offset.first);
2726 121955875 : libmesh_assert_equal_to (vig, vg_and_offset.second);
2727 : #endif
2728 308607976 : const unsigned int n_comp = node.n_comp_group(sys_num,vg);
2729 :
2730 : // There is a potential problem with h refinement. Imagine a
2731 : // quad9 that has a linear FE on it. Then, on the hanging side,
2732 : // it can falsely identify a DOF at the mid-edge node. This is why
2733 : // we go through FEInterface instead of node.n_comp() directly.
2734 308607976 : const unsigned int nc =
2735 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2736 59009498 : is_inf ?
2737 82322 : FEInterface::n_dofs_at_node(fe_type, add_p_level*p_level, &elem, n) :
2738 : #endif
2739 309007166 : ndan (type, fe_type.order + add_p_level*p_level, n);
2740 :
2741 : // If this is a non-vertex on a hanging node with extra
2742 : // degrees of freedom, we use the non-vertex dofs (which
2743 : // come in reverse order starting from the end, to
2744 : // simplify p refinement)
2745 308607976 : if (extra_hanging_dofs && !elem.is_vertex(n))
2746 : {
2747 10577632 : const int dof_offset = n_comp - nc;
2748 :
2749 : // We should never have fewer dofs than necessary on a
2750 : // node unless we're getting indices on a parent element,
2751 : // and we should never need the indices on such a node
2752 10577632 : if (dof_offset < 0)
2753 : {
2754 0 : libmesh_assert(!elem.active());
2755 0 : di.resize(di.size() + nc, DofObject::invalid_id);
2756 : }
2757 : else
2758 26499945 : for (int i=int(n_comp)-1; i>=dof_offset; i--)
2759 : {
2760 6267748 : const dof_id_type d =
2761 15922313 : node.dof_number(sys_num, vg, vig, i, n_comp);
2762 6267748 : libmesh_assert_not_equal_to (d, DofObject::invalid_id);
2763 15922313 : field_dofs_functor(elem, n, v, di, d);
2764 : }
2765 : }
2766 : // If this is a vertex or an element without extra hanging
2767 : // dofs, our dofs come in forward order coming from the
2768 : // beginning
2769 : else
2770 : {
2771 : // We have a good component index only if it's being
2772 : // used on this FE type (nc) *and* it's available on
2773 : // this DofObject (n_comp).
2774 298030344 : const unsigned int good_nc = std::min(n_comp, nc);
2775 568391108 : for (unsigned int i=0; i!=good_nc; ++i)
2776 : {
2777 108945169 : const dof_id_type d =
2778 270135280 : node.dof_number(sys_num, vg, vig, i, n_comp);
2779 108945169 : libmesh_assert_not_equal_to (d, DofObject::invalid_id);
2780 108945169 : libmesh_assert_less (d, this->n_dofs());
2781 270360764 : field_dofs_functor(elem, n, v, di, d);
2782 : }
2783 :
2784 : // With fewer good component indices than we need, e.g.
2785 : // due to subdomain expansion, the remaining expected
2786 : // indices are marked invalid.
2787 298030344 : if (n_comp < nc)
2788 0 : for (unsigned int i=n_comp; i!=nc; ++i)
2789 0 : di.push_back(DofObject::invalid_id);
2790 : }
2791 : }
2792 :
2793 : // If there are any element-based DOF numbers, get them
2794 58238830 : const unsigned int nc = FEInterface::n_dofs_per_elem(fe_type, add_p_level*p_level, &elem);
2795 :
2796 : // We should never have fewer dofs than necessary on an
2797 : // element unless we're getting indices on a parent element
2798 : // (and we should never need those indices) or off-domain for a
2799 : // subdomain-restricted variable (where invalid_id is the
2800 : // correct thing to return)
2801 58238830 : if (nc != 0)
2802 : {
2803 1753750 : const unsigned int n_comp = elem.n_comp_group(sys_num,vg);
2804 4590609 : if (elem.n_systems() > sys_num && nc <= n_comp)
2805 : {
2806 26945900 : for (unsigned int i=0; i<nc; i++)
2807 : {
2808 8088805 : const dof_id_type d =
2809 22190997 : elem.dof_number(sys_num, vg, vig, i, n_comp);
2810 8088805 : libmesh_assert_not_equal_to (d, DofObject::invalid_id);
2811 :
2812 22355291 : field_dofs_functor(elem, invalid_uint, v, di, d);
2813 : }
2814 : }
2815 : else
2816 : {
2817 0 : libmesh_assert(!elem.active() || fe_type.family == LAGRANGE || fe_type.family == SUBDIVISION);
2818 0 : di.resize(di.size() + nc, DofObject::invalid_id);
2819 : }
2820 : }
2821 : }
2822 58265858 : }
2823 :
2824 :
2825 :
2826 : template <typename ScalarDofsFunctor, typename FieldDofsFunctor>
2827 26175243 : void DofMap::dof_indices (const Elem * const elem,
2828 : std::vector<dof_id_type> & di,
2829 : const unsigned int vn,
2830 : ScalarDofsFunctor scalar_dofs_functor,
2831 : FieldDofsFunctor field_dofs_functor,
2832 : int p_level) const
2833 : {
2834 : // We now allow elem==nullptr to request just SCALAR dofs
2835 : // libmesh_assert(elem);
2836 :
2837 : // dof_indices() is a relatively light-weight function that is
2838 : // called millions of times in normal codes. Therefore, it is not a
2839 : // good candidate for logging, since the cost of the logging code
2840 : // itself is roughly on par with the time required to call
2841 : // dof_indices().
2842 : // LOG_SCOPE("dof_indices()", "DofMap");
2843 :
2844 : // Clear the DOF indices vector
2845 11015728 : di.clear();
2846 :
2847 : // Use the default p refinement level?
2848 26175243 : if (p_level == -12345)
2849 25763569 : p_level = elem ? elem->p_level() : 0;
2850 :
2851 26238597 : const unsigned int vg = this->_variable_group_numbers[vn];
2852 11015728 : const VariableGroup & var = this->variable_group(vg);
2853 26175243 : const unsigned int vig = vn - var.number();
2854 :
2855 : #ifdef DEBUG
2856 : // Check that sizes match in DEBUG mode
2857 11015728 : std::size_t tot_size = 0;
2858 : #endif
2859 :
2860 26175243 : if (elem && elem->type() == TRI3SUBDIVISION)
2861 : {
2862 : // Subdivision surface FE require the 1-ring around elem
2863 2712 : const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
2864 :
2865 : // Ghost subdivision elements have no real dofs
2866 8136 : if (!sd_elem->is_ghost())
2867 : {
2868 : // Determine the nodes contributing to element elem
2869 4608 : std::vector<const Node *> elem_nodes;
2870 6912 : MeshTools::Subdivision::find_one_ring(sd_elem, elem_nodes);
2871 :
2872 6912 : _dof_indices(*elem, p_level, di, vg, vig, elem_nodes.data(),
2873 : cast_int<unsigned int>(elem_nodes.size()), vn,
2874 : #ifdef DEBUG
2875 : tot_size,
2876 : #endif
2877 : field_dofs_functor);
2878 : }
2879 :
2880 8136 : return;
2881 : }
2882 :
2883 : // Get the dof numbers
2884 26265062 : if (var.type().family == SCALAR &&
2885 280285 : (!elem ||
2886 280317 : var.active_on_subdomain(elem->subdomain_id())))
2887 : {
2888 : #ifdef DEBUG
2889 97955 : tot_size += var.type().order;
2890 : #endif
2891 195910 : std::vector<dof_id_type> di_new;
2892 280317 : this->SCALAR_dof_indices(di_new,vn);
2893 280317 : scalar_dofs_functor(*elem, di, di_new);
2894 : }
2895 25886790 : else if (elem)
2896 25886790 : _dof_indices(*elem, p_level, di, vg, vig, elem->get_nodes(),
2897 25886790 : elem->n_nodes(), vn,
2898 : #ifdef DEBUG
2899 : tot_size,
2900 : #endif
2901 : field_dofs_functor);
2902 :
2903 : #ifdef DEBUG
2904 11013016 : libmesh_assert_equal_to (tot_size, di.size());
2905 : #endif
2906 : }
2907 :
2908 : inline
2909 18 : StaticCondensationDofMap & DofMap::get_static_condensation()
2910 : {
2911 18 : libmesh_assert(_sc);
2912 18 : return *_sc;
2913 : }
2914 :
2915 : inline
2916 28 : const StaticCondensationDofMap & DofMap::get_static_condensation() const
2917 : {
2918 28 : libmesh_assert(_sc);
2919 28 : return *_sc;
2920 : }
2921 :
2922 : inline
2923 0 : const ConstraintEliminationDofMap & DofMap::get_constraint_elimination() const
2924 : {
2925 0 : libmesh_assert(_constraint_elimination);
2926 0 : return *_constraint_elimination;
2927 : }
2928 :
2929 : inline const std::pair<unsigned int, unsigned int> &
2930 40 : DofMap::get_variable_array(const unsigned int vi) const
2931 : {
2932 24 : auto it = std::upper_bound(
2933 : _array_variables.begin(),
2934 : _array_variables.end(),
2935 : vi,
2936 72 : [](unsigned int value, const std::pair<unsigned int, unsigned int> & b) { return value < b.first; });
2937 :
2938 16 : libmesh_assert_msg(it != _array_variables.begin(),
2939 : "Passed in " << std::to_string(vi) << " is not in any of our array variables");
2940 16 : --it;
2941 16 : libmesh_assert_msg(vi < it->second,
2942 : "Passed in " << std::to_string(vi) << " is not in any of our array variables");
2943 56 : return *it;
2944 : }
2945 :
2946 : template <typename DofIndicesFunctor>
2947 56 : void DofMap::array_dof_indices(const DofIndicesFunctor & functor,
2948 : std::vector<dof_id_type> & di,
2949 : const unsigned int vn) const
2950 : {
2951 72 : const auto [begin, end] = this->get_variable_array(vn);
2952 40 : functor(di, begin);
2953 :
2954 56 : const unsigned int count = end - begin;
2955 : // We make count, which could be >> ntest, the inner index in hopes of vectorization
2956 56 : if (count > 1)
2957 : {
2958 40 : const dof_id_type component_size = di.size();
2959 56 : di.resize(count * component_size);
2960 :
2961 424 : const auto pack_container = [&di,
2962 : component_size](const unsigned int j,
2963 : const std::vector<dof_id_type> & j_dof_indices,
2964 32 : const unsigned int stride) {
2965 32 : if (&j_dof_indices != &di)
2966 4 : libmesh_assert(j_dof_indices.size() == component_size);
2967 672 : for (const auto i : make_range(component_size))
2968 768 : di[j * component_size + i] = j_dof_indices[i] + stride * j;
2969 : };
2970 40 : pack_container(0, di, 0);
2971 :
2972 56 : const auto & fe_type = _variable_groups[libmesh_map_find(_var_to_vg, vn)].type();
2973 72 : if (const bool lagrange = fe_type.family == LAGRANGE;
2974 56 : lagrange || (FEInterface::get_continuity(fe_type) == DISCONTINUOUS))
2975 : {
2976 42 : const auto stride = lagrange ? 1 : component_size;
2977 84 : for (const auto j : make_range((unsigned int)1, count))
2978 30 : pack_container(j, di, stride);
2979 : }
2980 : else
2981 : {
2982 14 : static thread_local std::vector<dof_id_type> work_dof_indices;
2983 4 : unsigned int j = 1;
2984 28 : for (const auto i : make_range(begin + 1, end))
2985 : {
2986 10 : functor(work_dof_indices, i);
2987 14 : pack_container(j++, work_dof_indices, 0);
2988 : }
2989 : }
2990 : }
2991 56 : }
2992 :
2993 : inline
2994 10548619 : unsigned int DofMap::n_vars() const
2995 : {
2996 21075704 : return cast_int<unsigned int>(_variables.size());
2997 : }
2998 :
2999 : inline
3000 2693086 : const std::string & DofMap::variable_name (const unsigned int i) const
3001 : {
3002 2693086 : libmesh_assert_less (i, _variables.size());
3003 :
3004 10727156 : return _variables[i].name();
3005 : }
3006 :
3007 : inline
3008 1368 : bool DofMap::identify_variable_groups () const
3009 : {
3010 3268 : return _identify_variable_groups;
3011 : }
3012 :
3013 : inline
3014 0 : void DofMap::identify_variable_groups (const bool ivg)
3015 : {
3016 0 : _identify_variable_groups = ivg;
3017 0 : }
3018 :
3019 : inline
3020 26092 : unsigned int DofMap::n_components(const MeshBase & mesh) const
3021 : {
3022 36604 : if (_variables.empty())
3023 10136 : return 0;
3024 :
3025 376 : const Variable & last = _variables.back();
3026 1264 : return last.first_scalar_number() + last.n_components(mesh);
3027 : }
3028 :
3029 : inline
3030 : unsigned int
3031 541077 : DofMap::variable_scalar_number (unsigned int var_num,
3032 : unsigned int component) const
3033 : {
3034 2511418 : return _variables[var_num].first_scalar_number() + component;
3035 : }
3036 :
3037 : inline
3038 3242 : const FEType & DofMap::variable_type (std::string_view var) const
3039 : {
3040 4152 : return _variables[this->variable_number(var)].type();
3041 : }
3042 :
3043 898 : inline bool DofMap::has_variable(std::string_view var) const
3044 : {
3045 898 : return _variable_numbers.count(var);
3046 : }
3047 :
3048 374373 : inline unsigned int DofMap::variable_number(std::string_view var) const
3049 : {
3050 1325387 : auto var_num = libmesh_map_find(_variable_numbers, var);
3051 374373 : libmesh_assert_equal_to(_variables[var_num].name(), var);
3052 374373 : return var_num;
3053 : }
3054 :
3055 : } // namespace libMesh
3056 :
3057 : #endif // LIBMESH_DOF_MAP_H
|
