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