Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : #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 230717 : 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 230581 : 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 230581 : 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 6774 : 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 requirements by matrices
274 : */
275 : void full_sparsity_pattern_needed();
276 :
277 : /**
278 : * Returns true iff the current policy when constructing sparsity
279 : * patterns is to explicitly account for sparsity entries created by
280 : * constraint matrix pre- and post- application.
281 : */
282 : bool constrained_sparsity_construction();
283 :
284 : /**
285 : * Clears the sparsity pattern
286 : */
287 : void clear_sparsity();
288 :
289 : /**
290 : * Remove any default ghosting functor(s). User-added ghosting
291 : * functors will be unaffected.
292 : *
293 : * Unless user-added equivalent ghosting functors exist, removing
294 : * the default coupling functor is only safe for explicit solves,
295 : * and removing the default algebraic ghosting functor is only safe
296 : * for codes where no evaluations on neighbor cells (e.g. no jump
297 : * error estimators) are done.
298 : *
299 : * Defaults can be restored manually via add_default_ghosting(), or
300 : * automatically if clear() returns the DofMap to a default state.
301 : */
302 : void remove_default_ghosting();
303 :
304 : /**
305 : * Add the default functor(s) for coupling and algebraic ghosting.
306 : * User-added ghosting functors will be unaffected.
307 : */
308 : void add_default_ghosting();
309 :
310 : /**
311 : * Iterator type for coupling and algebraic ghosting functor ranges.
312 : * This has changed in the past and may change again; code should
313 : * use auto or the type here.
314 : */
315 : typedef std::vector<GhostingFunctor *>::const_iterator GhostingFunctorIterator;
316 :
317 : /**
318 : * Adds a functor which can specify coupling requirements for
319 : * creation of sparse matrices.
320 : * Degree of freedom pairs which match the elements and variables
321 : * returned by these functors will be added to the sparsity pattern,
322 : * and the degrees of freedom which live on other processors will be
323 : * added to the send_list for use on ghosted vectors, and the
324 : * elements which live on other processors will be ghosted on a
325 : * distributed mesh.
326 : *
327 : * GhostingFunctor memory must be managed by the code which calls
328 : * this function; the GhostingFunctor lifetime is expected to extend
329 : * until either the functor is removed or the DofMap is destructed.
330 : *
331 : * When \p to_mesh is true, the \p coupling_functor is also added to
332 : * our associated mesh, to ensure that coupled elements do not get
333 : * lost during mesh distribution. (if coupled elements were
334 : * *already* lost there's no getting them back after the fact,
335 : * sorry)
336 : *
337 : * If \p to_mesh is false, no change to mesh ghosting is made;
338 : * the Mesh must already have ghosting functor(s) specifying a
339 : * superset of \p coupling_functor or this is a horrible bug.
340 : */
341 : void add_coupling_functor(GhostingFunctor & coupling_functor,
342 : bool to_mesh = true);
343 :
344 : /**
345 : * Adds a functor which can specify coupling requirements for
346 : * creation of sparse matrices.
347 : *
348 : * GhostingFunctor memory when using this method is managed by the
349 : * shared_ptr mechanism.
350 : */
351 : void add_coupling_functor(std::shared_ptr<GhostingFunctor> coupling_functor,
352 : bool to_mesh = true)
353 : { _shared_functors[coupling_functor.get()] = coupling_functor;
354 : this->add_coupling_functor(*coupling_functor, to_mesh); }
355 :
356 : /**
357 : * Removes a functor which was previously added to the set of
358 : * coupling functors, from both this DofMap and from the underlying
359 : * mesh.
360 : */
361 : void remove_coupling_functor(GhostingFunctor & coupling_functor);
362 :
363 : /**
364 : * Beginning of range of coupling functors
365 : */
366 314002 : GhostingFunctorIterator coupling_functors_begin() const
367 637556 : { return _coupling_functors.begin(); }
368 :
369 : /**
370 : * End of range of coupling functors
371 : */
372 314002 : GhostingFunctorIterator coupling_functors_end() const
373 637556 : { return _coupling_functors.end(); }
374 :
375 : /**
376 : * Default coupling functor
377 : */
378 0 : DefaultCoupling & default_coupling() { return *_default_coupling; }
379 :
380 : /**
381 : * Adds a functor which can specify algebraic ghosting requirements
382 : * for use with distributed vectors. Degrees of freedom on other
383 : * processors which match the elements and variables returned by
384 : * these functors will be added to the send_list, and the elements
385 : * on other processors will be ghosted on a distributed mesh, so
386 : * that the elements can always be found and the solutions on them
387 : * will always be evaluable.
388 : *
389 : * GhostingFunctor memory must be managed by the code which calls
390 : * this function; the GhostingFunctor lifetime is expected to extend
391 : * until either the functor is removed or the DofMap is destructed.
392 : *
393 : * When \p to_mesh is true, the \p coupling_functor is also added to
394 : * our associated mesh, to ensure that evaluable elements do not get
395 : * lost during mesh distribution. (if evaluable elements were
396 : * *already* lost there's no getting them back after the fact,
397 : * sorry)
398 : *
399 : * If \p to_mesh is false, no change to mesh ghosting is made;
400 : * the Mesh must already have ghosting functor(s) specifying a
401 : * superset of \p evaluable_functor or this is a horrible bug.
402 : */
403 : void add_algebraic_ghosting_functor(GhostingFunctor & evaluable_functor,
404 : bool to_mesh = true);
405 :
406 : /**
407 : * Adds a functor which can specify algebraic ghosting requirements
408 : * for use with distributed vectors.
409 : *
410 : * GhostingFunctor memory when using this method is managed by the
411 : * shared_ptr mechanism.
412 : */
413 : void add_algebraic_ghosting_functor(std::shared_ptr<GhostingFunctor> evaluable_functor,
414 : bool to_mesh = true)
415 : { _shared_functors[evaluable_functor.get()] = evaluable_functor;
416 : this->add_algebraic_ghosting_functor(*evaluable_functor, to_mesh); }
417 :
418 : /**
419 : * Removes a functor which was previously added to the set of
420 : * algebraic ghosting functors, from both this DofMap and from the
421 : * underlying mesh.
422 : */
423 : void remove_algebraic_ghosting_functor(GhostingFunctor & evaluable_functor);
424 :
425 : /**
426 : * Beginning of range of algebraic ghosting functors
427 : */
428 1568 : GhostingFunctorIterator algebraic_ghosting_functors_begin() const
429 9180 : { return _algebraic_ghosting_functors.begin(); }
430 :
431 : /**
432 : * End of range of algebraic ghosting functors
433 : */
434 1568 : GhostingFunctorIterator algebraic_ghosting_functors_end() const
435 9180 : { return _algebraic_ghosting_functors.end(); }
436 :
437 : /**
438 : * Default algebraic ghosting functor
439 : */
440 0 : DefaultCoupling & default_algebraic_ghosting() { return *_default_evaluating; }
441 :
442 : /**
443 : * Attach an object to use to populate the
444 : * sparsity pattern with extra entries.
445 : *
446 : * Care must be taken that when adding entries they are sorted into the Rows
447 : *
448 : * Further, you _must_ modify n_nz and n_oz properly!
449 : *
450 : * This is an advanced function... use at your own peril!
451 : */
452 : void attach_extra_sparsity_object (SparsityPattern::AugmentSparsityPattern & asp)
453 : {
454 : _augment_sparsity_pattern = &asp;
455 : }
456 :
457 : /**
458 : * Attach a function pointer to use as a callback to populate the
459 : * sparsity pattern with extra entries.
460 : *
461 : * Care must be taken that when adding entries they are sorted into the Rows
462 : *
463 : * Further, you _must_ modify n_nz and n_oz properly!
464 : *
465 : * This is an advanced function... use at your own peril!
466 : */
467 : void attach_extra_sparsity_function(void (*func)(SparsityPattern::Graph & sparsity,
468 : std::vector<dof_id_type> & n_nz,
469 : std::vector<dof_id_type> & n_oz,
470 : void *),
471 : void * context = nullptr)
472 : { _extra_sparsity_function = func; _extra_sparsity_context = context; }
473 :
474 : /**
475 : * Attach an object to populate the send_list with extra entries.
476 : * This should only add to the send list, but no checking is done
477 : * to enforce this behavior.
478 : *
479 : * This is an advanced function... use at your own peril!
480 : */
481 : void attach_extra_send_list_object (DofMap::AugmentSendList & asl)
482 : {
483 : _augment_send_list = &asl;
484 : }
485 :
486 : /**
487 : * Attach a function pointer to use as a callback to populate the
488 : * send_list with extra entries.
489 : */
490 : void attach_extra_send_list_function(void (*func)(std::vector<dof_id_type> &, void *),
491 : void * context = nullptr)
492 : { _extra_send_list_function = func; _extra_send_list_context = context; }
493 :
494 : /**
495 : * Takes the \p _send_list vector (which may have duplicate entries)
496 : * and sorts it. The duplicate entries are then removed, resulting in
497 : * a sorted \p _send_list with unique entries. Also calls any user-provided
498 : * methods for adding to the send list.
499 : */
500 : void prepare_send_list ();
501 :
502 : /**
503 : * Clears the \p _send_list vector. This should be done in order to completely
504 : * rebuild the send_list from scratch rather than merely adding to the existing
505 : * send_list.
506 : */
507 14412 : void clear_send_list ()
508 : {
509 14412 : _send_list.clear();
510 14412 : }
511 :
512 : /**
513 : * Clears the \p _send_list vector and then rebuilds it. This may be needed
514 : * in special situations, for example when an algebraic coupling functor cannot
515 : * be added to the \p DofMap until after it is completely setup. Then this method
516 : * can be used to rebuild the send_list once the algebraic coupling functor is
517 : * added. Note that while this will recommunicate constraints with the updated
518 : * send_list, this does assume no new constraints have been added since the previous
519 : * reinit_constraints call.
520 : */
521 : void reinit_send_list (MeshBase & mesh);
522 :
523 :
524 : /**
525 : * \returns A constant reference to the \p _send_list for this processor.
526 : *
527 : * The \p _send_list contains the global indices of all the
528 : * variables in the global solution vector that influence the
529 : * current processor. This information can be used for gathers at
530 : * each solution step to retrieve solution values needed for
531 : * computation.
532 : */
533 2601597 : const std::vector<dof_id_type> & get_send_list() const { return _send_list; }
534 :
535 : /**
536 : * \returns A constant reference to the \p _n_nz list for this processor.
537 : *
538 : * The vector contains the bandwidth of the on-processor coupling for each
539 : * row of the global matrix that the current processor owns. This
540 : * information can be used to preallocate space for a parallel sparse matrix.
541 : */
542 : const std::vector<dof_id_type> & get_n_nz() const
543 : {
544 : libmesh_assert(_sp);
545 : return _sp->get_n_nz();
546 : }
547 :
548 : /**
549 : * \returns A constant reference to the \p _n_oz list for this processor.
550 : *
551 : * The vector contains the bandwidth of the off-processor coupling for each
552 : * row of the global matrix that the current processor owns. This
553 : * information can be used to preallocate space for a parallel sparse matrix.
554 : */
555 : const std::vector<dof_id_type> & get_n_oz() const
556 : {
557 : libmesh_assert(_sp);
558 : return _sp->get_n_oz();
559 : }
560 :
561 :
562 : /**
563 : * \returns A constant pointer to the sparsity pattern stored here,
564 : * once that has been computed. Returns null if no sparsity pattern
565 : * has yet been computed.
566 : *
567 : * If need_full_sparsity_pattern is false, the "sparsity pattern"
568 : * may only own n_nz and n_oz lists.
569 : */
570 638 : const SparsityPattern::Build * get_sparsity_pattern() const
571 : {
572 638 : return _sp.get();
573 : }
574 :
575 : // /**
576 : // * Add an unknown of order \p order and finite element type
577 : // * \p type to the system of equations.
578 : // */
579 : // void add_variable (const Variable & var);
580 :
581 : /**
582 : * Add a group of unknowns of order \p order and finite element type
583 : * \p type to the system of equations.
584 : */
585 : void add_variable_group (VariableGroup var_group);
586 :
587 : /**
588 : * Specify whether or not we perform an extra (opt-mode enabled) check
589 : * for constraint loops. If a constraint loop is present then
590 : * the system constraints are not valid, so if \p error_on_constraint_loop
591 : * is true we will throw an error in this case.
592 : *
593 : * \note We previously referred to these types of constraints as
594 : * "cyclic" but that has now been deprecated, and these will now
595 : * instead be referred to as "constraint loops" in libMesh.
596 : */
597 : void set_error_on_cyclic_constraint(bool error_on_cyclic_constraint);
598 : void set_error_on_constraint_loop(bool error_on_constraint_loop);
599 :
600 : /**
601 : * \returns The \p VariableGroup description object for group \p g.
602 : */
603 : const VariableGroup & variable_group (const unsigned int c) const;
604 :
605 : const Variable & variable (const unsigned int c) const override;
606 :
607 : /**
608 : * \returns The approximation order for variable \p c.
609 : */
610 : Order variable_order (const unsigned int c) const;
611 :
612 : /**
613 : * \returns The approximation order for \p VariableGroup \p vg.
614 : */
615 : Order variable_group_order (const unsigned int vg) const;
616 :
617 : /**
618 : * \returns The finite element type for variable \p c.
619 : */
620 : const FEType & variable_type (const unsigned int c) const;
621 :
622 : /**
623 : * \returns The finite element type for \p VariableGroup \p vg.
624 : */
625 : const FEType & variable_group_type (const unsigned int vg) const;
626 :
627 : /**
628 : * \returns The number of variables in the global solution vector. Defaults
629 : * to 1, should be 1 for a scalar equation, 3 for 2D incompressible Navier
630 : * Stokes (u,v,p), etc...
631 : */
632 11769560 : unsigned int n_variable_groups() const
633 23330258 : { return cast_int<unsigned int>(_variable_groups.size()); }
634 :
635 3205066 : unsigned int n_variables() const override
636 3271440 : { return cast_int<unsigned int>(_variables.size()); }
637 :
638 : /**
639 : * \returns The variable group number that the provided variable number belongs to
640 : */
641 : unsigned int var_group_from_var_number(unsigned int var_num) const;
642 :
643 : /**
644 : * \returns \p true if the variables are capable of being stored in a blocked
645 : * form. Presently, this means that there can only be one variable group,
646 : * and that the group has more than one variable.
647 : */
648 39892 : bool has_blocked_representation() const
649 : {
650 41206 : return ((this->n_variable_groups() == 1) && (this->n_variables() > 1));
651 : }
652 :
653 : /**
654 : * \returns The block size, if the variables are amenable to block storage.
655 : * Otherwise 1.
656 : * This routine was originally designed to enable a blocked storage, but
657 : * it turns out this information is still super useful for solvers even when
658 : * we do not use the blocked storage (e.g., MATMPIBAIJ in PETSc). For example (in PCHMG),
659 : * for a system of PDEs, to construct an efficient multilevel preconditioner, we coarsen
660 : * the matrix of one single PDE instead of the entire huge matrix. In order to
661 : * accomplish this, we need to know many PDEs we have. Another use case,
662 : * the fieldsplit preconditioner can be constructed in place with this info without
663 : * involving any user efforts.
664 : */
665 41206 : unsigned int block_size() const
666 : {
667 39982 : return (this->has_blocked_representation() ? this->n_variables() : 1);
668 : }
669 :
670 : using DofMapBase::n_dofs;
671 : /**
672 : * \returns The total number of degrees of freedom for a particular
673 : * variable \p vn.
674 : */
675 : dof_id_type n_dofs(const unsigned int vn) const
676 : {
677 : dof_id_type n = this->n_local_dofs(vn);
678 : this->comm().sum(n);
679 : return n;
680 : }
681 :
682 : /**
683 : * \returns The number of SCALAR dofs.
684 : */
685 255348 : dof_id_type n_SCALAR_dofs() const { return _n_SCALAR_dofs; }
686 :
687 : using DofMapBase::n_local_dofs;
688 : /**
689 : * \returns The number of degrees of freedom on this processor for a
690 : * particular variable \p vn. This is an O(N) operation on serialized or
691 : * O(N/Nproc) operation on distributed meshes.
692 : */
693 : dof_id_type n_local_dofs(const unsigned int vn) const
694 : {
695 : dof_id_type n;
696 : this->local_variable_indices(n, _mesh, vn);
697 : return n;
698 : }
699 :
700 : /**
701 : * \returns The number of degrees of freedom on each partition for a
702 : * particular variable \p vn.
703 : */
704 : std::vector<dof_id_type> n_dofs_per_processor(const unsigned int vn) const
705 : {
706 : std::vector<dof_id_type> n_local_dofs(this->n_processors(), 0);
707 : this->comm().allgather(this->n_local_dofs(vn), n_local_dofs);
708 : return n_local_dofs;
709 : }
710 :
711 : /**
712 : * \returns The processor id that owns the dof index \p dof
713 : */
714 94065 : processor_id_type dof_owner(const dof_id_type dof) const
715 : { std::vector<dof_id_type>::const_iterator ub =
716 94065 : std::upper_bound(_end_df.begin(), _end_df.end(), dof);
717 2602 : libmesh_assert (ub != _end_df.end());
718 94065 : return cast_int<processor_id_type>(ub - _end_df.begin());
719 : }
720 :
721 : void dof_indices (const Elem * const elem,
722 : std::vector<dof_id_type> & di) const;
723 :
724 : /**
725 : * Fills the vector \p di with the global degree of freedom indices
726 : * for the element. For one variable, and potentially for a
727 : * non-default element p refinement level
728 : */
729 : void dof_indices (const Elem * const elem,
730 : std::vector<dof_id_type> & di,
731 : const unsigned int vn,
732 : int p_level = -12345) const override;
733 :
734 : /**
735 : * Retrieves degree of freedom indices for a given \p elem and then performs actions for these
736 : * indices defined by the user-provided functors \p scalar_dofs_functor and \p field_dofs_functor.
737 : * This API is useful when a user wants to do more than simply fill a degree of freedom container
738 : * @param elem The element to get degrees of freedom for
739 : * @param di A container for degrees of freedom. It is up to the provided functors how this gets
740 : * filled
741 : * @param vn The variable number to retrieve degrees of freedom for
742 : * @param scalar_dofs_functor The functor that acts on scalar degrees of freedom. This functor has
743 : * the interface:
744 : * void scalar_dofs_functor(const Elem & elem,
745 : * std::vector<dof_id_type> & di,
746 : * const std::vector<dof_id_type> & scalar_dof_indices)
747 : * where \p di is the degree of freedom container described above and
748 : * \p scalar_dof_indices are the scalar dof indices available to
749 : * \p elem
750 : * @param field_dofs_functor The functor that acts on "field" (e.g. non-scalar, non-global)
751 : * degrees of freedom. This functor has
752 : * the interface:
753 : * void field_dofs_functor(const Elem & elem,
754 : * const unsigned int node_num,
755 : * const unsigned int var_num,
756 : * std::vector<dof_id_type> & di,
757 : * const dof_id_type field_dof)
758 : * where \p field_dof represents a field degree of freedom to act on and
759 : * is associated with \p node_num and \p var_num. If the degree of
760 : * freedom is elemental than \p node_num will be \p invalid_uint. \p di
761 : * is again the degree of freedom container provided above
762 : */
763 : template <typename ScalarDofsFunctor, typename FieldDofsFunctor>
764 : void dof_indices(const Elem * const elem,
765 : std::vector<dof_id_type> & di,
766 : const unsigned int vn,
767 : ScalarDofsFunctor scalar_dofs_functor,
768 : FieldDofsFunctor field_dofs_functor,
769 : int p_level = -12345) const;
770 :
771 : /**
772 : * Fills the vector \p di with the global degree of freedom indices
773 : * for the \p node.
774 : */
775 : void dof_indices (const Node * const node,
776 : std::vector<dof_id_type> & di) const;
777 :
778 : /**
779 : * Fills the vector \p di with the global degree of freedom indices
780 : * for the \p node, for one variable \p vn.
781 : */
782 : void dof_indices (const Node * const node,
783 : std::vector<dof_id_type> & di,
784 : const unsigned int vn) const override;
785 :
786 : /**
787 : * Appends to the vector \p di the global degree of freedom indices
788 : * for \p elem.node_ref(n), for one variable \p vn. On hanging
789 : * nodes with both vertex and non-vertex DoFs, only those indices
790 : * which are directly supported on \p elem are included.
791 : */
792 : void dof_indices (const Elem & elem,
793 : unsigned int n,
794 : std::vector<dof_id_type> & di,
795 : const unsigned int vn) const;
796 :
797 : #ifdef LIBMESH_ENABLE_AMR
798 :
799 : /**
800 : * Appends to the vector \p di the old global degree of freedom
801 : * indices for \p elem.node_ref(n), for one variable \p vn. On
802 : * hanging nodes with both vertex and non-vertex DoFs, only those
803 : * indices which are directly supported on \p elem are included.
804 : */
805 : void old_dof_indices (const Elem & elem,
806 : unsigned int n,
807 : std::vector<dof_id_type> & di,
808 : const unsigned int vn) const;
809 :
810 : #endif // LIBMESH_ENABLE_AMR
811 :
812 : /**
813 : * Fills the vector \p di with the global degree of freedom indices
814 : * corresponding to the SCALAR variable vn. If old_dofs=true,
815 : * the old SCALAR dof indices are returned.
816 : *
817 : * \note We do not need to pass in an element since SCALARs are
818 : * global variables.
819 : */
820 : void SCALAR_dof_indices (std::vector<dof_id_type> & di,
821 : const unsigned int vn,
822 : const bool old_dofs=false) const;
823 :
824 : /**
825 : * \returns \p true if degree of freedom index \p dof_index
826 : * is either a local index or in the \p send_list.
827 : *
828 : * \note This is an O(logN) operation for a send_list of size N; we
829 : * don't cache enough information for O(1) right now.
830 : */
831 : bool semilocal_index (dof_id_type dof_index) const;
832 :
833 : /**
834 : * \returns \p true if all degree of freedom indices in \p
835 : * dof_indices are either local indices or in the \p send_list.
836 : *
837 : * \note This is an O(logN) operation for a send_list of size N; we
838 : * don't cache enough information for O(1) right now.
839 : */
840 : bool all_semilocal_indices (const std::vector<dof_id_type> & dof_indices) const;
841 :
842 : /**
843 : * \returns \p true if degree of freedom index \p dof_index
844 : * is a local index.
845 : */
846 72258077 : bool local_index (dof_id_type dof_index) const
847 78434465 : { return (dof_index >= this->first_dof()) && (dof_index < this->end_dof()); }
848 :
849 : /**
850 : * \returns \p true iff our solutions can be locally evaluated on
851 : * \p obj (which should be an Elem or a Node) for variable number \p
852 : * var_num (for all variables, if \p var_num is invalid_uint)
853 : */
854 : template <typename DofObjectSubclass>
855 : bool is_evaluable(const DofObjectSubclass & obj,
856 : unsigned int var_num = libMesh::invalid_uint) const;
857 :
858 : /**
859 : * Allow the implicit_neighbor_dofs flag to be set programmatically.
860 : * This overrides the --implicit_neighbor_dofs commandline option.
861 : * We can use this to set the implicit neighbor dofs option differently
862 : * for different systems, whereas the commandline option is the same
863 : * for all systems.
864 : */
865 : void set_implicit_neighbor_dofs(bool implicit_neighbor_dofs);
866 :
867 : /**
868 : * Set the _verify_dirichlet_bc_consistency flag.
869 : */
870 : void set_verify_dirichlet_bc_consistency(bool val);
871 :
872 : /**
873 : * Tells other library functions whether or not this problem
874 : * includes coupling between dofs in neighboring cells, as can
875 : * currently be specified on the command line or inferred from
876 : * the use of all discontinuous variables.
877 : */
878 : bool use_coupled_neighbor_dofs(const MeshBase & mesh) const;
879 :
880 : /**
881 : * Builds the local element vector \p Ue from the global vector \p Ug,
882 : * accounting for any constrained degrees of freedom. For an element
883 : * without constrained degrees of freedom this is the trivial mapping
884 : * \f$ Ue[i] = Ug[dof_indices[i]] \f$
885 : *
886 : * \note The user must ensure that the element vector \p Ue is
887 : * properly sized when calling this method. This is because there
888 : * is no \p resize() method in the \p DenseVectorBase<> class.
889 : */
890 : void extract_local_vector (const NumericVector<Number> & Ug,
891 : const std::vector<dof_id_type> & dof_indices,
892 : DenseVectorBase<Number> & Ue) const;
893 :
894 : /**
895 : * If T == dof_id_type, counts, if T == std::vector<dof_id_type>, fills an
896 : * array of, those dof indices which belong to the given variable number and
897 : * live on the current processor.
898 : */
899 : template <typename T, std::enable_if_t<std::is_same_v<T, dof_id_type> ||
900 : std::is_same_v<T, std::vector<dof_id_type>>, int> = 0>
901 : void local_variable_indices(T & idx,
902 : const MeshBase & mesh,
903 : unsigned int var_num) const;
904 :
905 : /**
906 : * If T == dof_id_type, counts, if T == std::vector<dof_id_type>, fills an
907 : * array of, those dof indices which belong to the given variable number and
908 : * live on the current processor.
909 : */
910 : template <typename T,
911 : std::enable_if_t<std::is_same_v<T, dof_id_type> ||
912 : std::is_same_v<T, std::vector<dof_id_type>>,
913 : int> = 0>
914 : void local_variable_indices(T & idx, unsigned int var_num) const
915 : { this->local_variable_indices(idx, this->_mesh, var_num); }
916 :
917 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
918 :
919 : //--------------------------------------------------------------------
920 : // Constraint-specific methods
921 : /**
922 : * \returns The total number of constrained degrees of freedom
923 : * in the problem.
924 : */
925 : dof_id_type n_constrained_dofs() const;
926 :
927 : /**
928 : * \returns The number of constrained degrees of freedom
929 : * on this processor.
930 : */
931 : dof_id_type n_local_constrained_dofs() const;
932 :
933 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
934 : /**
935 : * \returns The total number of constrained Nodes
936 : * in the mesh.
937 : */
938 : dof_id_type n_constrained_nodes() const
939 : { return cast_int<dof_id_type>(_node_constraints.size()); }
940 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
941 :
942 : /**
943 : * Rebuilds the raw degree of freedom and DofObject constraints,
944 : * based on attached DirichletBoundary objects and on non-conforming
945 : * interface in adapted meshes.
946 : *
947 : * A time is specified for use in building time-dependent Dirichlet
948 : * constraints.
949 : */
950 : void create_dof_constraints (const MeshBase &, Real time=0);
951 :
952 : /**
953 : * Gathers constraint equation dependencies from other processors
954 : */
955 : void allgather_recursive_constraints (MeshBase &);
956 :
957 : /**
958 : * Sends constraint equations to constraining processors
959 : */
960 : void scatter_constraints (MeshBase &);
961 :
962 : /**
963 : * Helper function for querying about constraint equations on other
964 : * processors. If any id in \p requested_dof_ids is constrained on
965 : * another processor, its constraint will be added on this processor
966 : * as well. If \p look_for_constrainees is true, then constraints
967 : * will also be returned if the id appears as a constraining value
968 : * not just if it appears as a constrained value.
969 : *
970 : * This function operates recursively: if the constraint for a
971 : * constrained dof is newly added locally, then any other dofs which
972 : * constrain it are queried to see if they are in turn constrained,
973 : * and so on.
974 : */
975 : void gather_constraints (MeshBase & mesh,
976 : std::set<dof_id_type> & unexpanded_dofs,
977 : bool look_for_constrainees);
978 :
979 : /**
980 : * Postprocesses any constrained degrees of freedom
981 : * to be constrained only in terms of unconstrained dofs, then adds
982 : * unconstrained dofs to the send_list and prepares that for use.
983 : * This should be run after both system (create_dof_constraints) and
984 : * user constraints have all been added.
985 : */
986 : void process_constraints (MeshBase &);
987 :
988 : /**
989 : * Throw an error if we detect any constraint loops, i.e.
990 : * A -> B -> C -> A
991 : * that is, "dof A is constrained in terms of dof B which is
992 : * constrained in terms of dof C which is constrained in terms of
993 : * dof A", since these are not supported by libMesh and give
994 : * erroneous results if they are present.
995 : *
996 : * \note The original "cyclic constraint" terminology was
997 : * unfortunate since the word cyclic is used by some software to
998 : * indicate an actual type of rotational/angular constraint and not
999 : * (as here) a cyclic graph. The former nomenclature will eventually
1000 : * be deprecated in favor of "constraint loop".
1001 : */
1002 : void check_for_cyclic_constraints();
1003 : void check_for_constraint_loops();
1004 :
1005 : /**
1006 : * Adds a copy of the user-defined row to the constraint matrix, using
1007 : * an inhomogeneous right-hand-side for the constraint equation.
1008 : */
1009 : void add_constraint_row (const dof_id_type dof_number,
1010 : const DofConstraintRow & constraint_row,
1011 : const Number constraint_rhs,
1012 : const bool forbid_constraint_overwrite);
1013 :
1014 : /**
1015 : * Adds a copy of the user-defined row to the constraint matrix,
1016 : * using an inhomogeneous right-hand-side for the adjoint constraint
1017 : * equation.
1018 : *
1019 : * \p forbid_constraint_overwrite here only tests for overwriting
1020 : * the rhs. This method should only be used when an equivalent
1021 : * constraint (with a potentially different rhs) already exists for
1022 : * the primal problem.
1023 : */
1024 : void add_adjoint_constraint_row (const unsigned int qoi_index,
1025 : const dof_id_type dof_number,
1026 : const DofConstraintRow & constraint_row,
1027 : const Number constraint_rhs,
1028 : const bool forbid_constraint_overwrite);
1029 :
1030 : /**
1031 : * Adds a copy of the user-defined row to the constraint matrix, using
1032 : * a homogeneous right-hand-side for the constraint equation.
1033 : * By default, produces an error if the DOF was already constrained.
1034 : */
1035 19370 : void add_constraint_row (const dof_id_type dof_number,
1036 : const DofConstraintRow & constraint_row,
1037 : const bool forbid_constraint_overwrite = true)
1038 232736 : { add_constraint_row(dof_number, constraint_row, 0., forbid_constraint_overwrite); }
1039 :
1040 : /**
1041 : * \returns An iterator pointing to the first DoF constraint row.
1042 : */
1043 : DofConstraints::const_iterator constraint_rows_begin() const
1044 : { return _dof_constraints.begin(); }
1045 :
1046 : /**
1047 : * \returns An iterator pointing just past the last DoF constraint row.
1048 : */
1049 : DofConstraints::const_iterator constraint_rows_end() const
1050 : { return _dof_constraints.end(); }
1051 :
1052 : /**
1053 : * Provide a const accessor to the DofConstraints map. This allows the user
1054 : * to quickly search the data structure rather than just iterating over it.
1055 : */
1056 2812 : const DofConstraints & get_dof_constraints() const { return _dof_constraints; }
1057 :
1058 : void stash_dof_constraints()
1059 : {
1060 : libmesh_assert(_stashed_dof_constraints.empty());
1061 : _dof_constraints.swap(_stashed_dof_constraints);
1062 : }
1063 :
1064 : void unstash_dof_constraints()
1065 : {
1066 : libmesh_assert(_dof_constraints.empty());
1067 : _dof_constraints.swap(_stashed_dof_constraints);
1068 : }
1069 :
1070 : /**
1071 : * Similar to the stash/unstash_dof_constraints() API, but swaps
1072 : * _dof_constraints and _stashed_dof_constraints without asserting
1073 : * that the source or destination is empty first.
1074 : *
1075 : * \note There is an implicit assumption that swapping between sets
1076 : * of Constraints does not change the sparsity pattern or expand the
1077 : * send_list, since the only thing changed is the DofConstraints
1078 : * themselves. This is intended to work for swapping between
1079 : * DofConstraints A and B, where A is used to define the send_list,
1080 : * and B is a subset of A.
1081 : */
1082 : void swap_dof_constraints()
1083 : {
1084 : _dof_constraints.swap(_stashed_dof_constraints);
1085 : }
1086 :
1087 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1088 : /**
1089 : * \returns An iterator pointing to the first Node constraint row.
1090 : */
1091 : NodeConstraints::const_iterator node_constraint_rows_begin() const
1092 : { return _node_constraints.begin(); }
1093 :
1094 : /**
1095 : * \returns An iterator pointing just past the last Node constraint row.
1096 : */
1097 : NodeConstraints::const_iterator node_constraint_rows_end() const
1098 : { return _node_constraints.end(); }
1099 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
1100 :
1101 : /**
1102 : * \returns \p true if the degree of freedom dof is constrained,
1103 : * \p false otherwise.
1104 : */
1105 : bool is_constrained_dof (const dof_id_type dof) const;
1106 :
1107 : /**
1108 : * \returns \p true if the system has any heterogeneous constraints for
1109 : * adjoint solution \p qoi_num, \p false otherwise.
1110 : */
1111 : bool has_heterogeneous_adjoint_constraints (const unsigned int qoi_num) const;
1112 :
1113 : /**
1114 : * Backwards compatibility with misspelling.
1115 : */
1116 4796 : bool has_heterogenous_adjoint_constraints (const unsigned int qoi_num) const
1117 : {
1118 83110 : return this->has_heterogeneous_adjoint_constraints (qoi_num);
1119 : }
1120 :
1121 : /**
1122 : * \returns The heterogeneous constraint value if the degree of
1123 : * freedom \p dof has a heterogeneous constraint for adjoint solution
1124 : * \p qoi_num, zero otherwise.
1125 : */
1126 : Number has_heterogeneous_adjoint_constraint (const unsigned int qoi_num,
1127 : const dof_id_type dof) const;
1128 :
1129 : /**
1130 : * Backwards compatibility with misspelling.
1131 : */
1132 1094202 : Number has_heterogenous_adjoint_constraint (const unsigned int qoi_num,
1133 : const dof_id_type dof) const
1134 : {
1135 1415182 : return this->has_heterogeneous_adjoint_constraint (qoi_num, dof);
1136 : }
1137 :
1138 : /**
1139 : * \returns A reference to the set of right-hand-side values in
1140 : * primal constraint equations
1141 : */
1142 : DofConstraintValueMap & get_primal_constraint_values();
1143 :
1144 : /**
1145 : * \returns \p true if the Node is constrained,
1146 : * false otherwise.
1147 : */
1148 : bool is_constrained_node (const Node * node) const;
1149 :
1150 : /**
1151 : * Prints (from processor 0) all DoF and Node constraints. If \p
1152 : * print_nonlocal is true, then each constraint is printed once for
1153 : * each processor that knows about it, which may be useful for \p
1154 : * DistributedMesh debugging.
1155 : */
1156 : void print_dof_constraints(std::ostream & os=libMesh::out,
1157 : bool print_nonlocal=false) const;
1158 :
1159 : /**
1160 : * Gets a string reporting all DoF and Node constraints local to
1161 : * this processor. If \p print_nonlocal is true, then nonlocal
1162 : * constraints which are locally known are included.
1163 : */
1164 : std::string get_local_constraints(bool print_nonlocal=false) const;
1165 :
1166 :
1167 : /**
1168 : * Tests the constrained degrees of freedom on the numeric vector \p v, which
1169 : * represents a solution defined on the mesh, returning a pair whose first
1170 : * entry is the maximum absolute error on a constrained DoF and whose second
1171 : * entry is the maximum relative error. Useful for debugging purposes.
1172 : *
1173 : * If \p v == nullptr, the system solution vector is tested.
1174 : */
1175 : std::pair<Real, Real> max_constraint_error(const System & system,
1176 : NumericVector<Number> * v = nullptr) const;
1177 :
1178 : #endif // LIBMESH_ENABLE_CONSTRAINTS
1179 :
1180 : //--------------------------------------------------------------------
1181 : // Constraint-specific methods
1182 : // Some of these methods are enabled (but inlined away to nothing)
1183 : // when constraints are disabled at configure-time. This is to
1184 : // increase API compatibility of user code with different library
1185 : // builds.
1186 :
1187 : /**
1188 : * Constrains the element matrix. This method requires the
1189 : * element matrix to be square, in which case the elem_dofs
1190 : * correspond to the global DOF indices of both the rows and
1191 : * columns of the element matrix. For this case the rows
1192 : * and columns of the matrix necessarily correspond to variables
1193 : * of the same approximation order.
1194 : *
1195 : * If \p asymmetric_constraint_rows is set to true (as it is by
1196 : * default), constraint row equations will be reinforced in a way
1197 : * which breaks matrix symmetry but makes inexact linear solver
1198 : * solutions more likely to satisfy hanging node constraints.
1199 : */
1200 : void constrain_element_matrix (DenseMatrix<Number> & matrix,
1201 : std::vector<dof_id_type> & elem_dofs,
1202 : bool asymmetric_constraint_rows = true) const;
1203 :
1204 : /**
1205 : * Constrains the element matrix. This method allows the
1206 : * element matrix to be non-square, in which case the row_dofs
1207 : * and col_dofs may be of different size and correspond to
1208 : * variables approximated in different spaces.
1209 : */
1210 : void constrain_element_matrix (DenseMatrix<Number> & matrix,
1211 : std::vector<dof_id_type> & row_dofs,
1212 : std::vector<dof_id_type> & col_dofs,
1213 : bool asymmetric_constraint_rows = true) const;
1214 :
1215 : /**
1216 : * Constrains the element vector.
1217 : */
1218 : void constrain_element_vector (DenseVector<Number> & rhs,
1219 : std::vector<dof_id_type> & dofs,
1220 : bool asymmetric_constraint_rows = true) const;
1221 :
1222 : /**
1223 : * Constrains the element matrix and vector. This method requires
1224 : * the element matrix to be square, in which case the elem_dofs
1225 : * correspond to the global DOF indices of both the rows and
1226 : * columns of the element matrix. For this case the rows
1227 : * and columns of the matrix necessarily correspond to variables
1228 : * of the same approximation order.
1229 : */
1230 : void constrain_element_matrix_and_vector (DenseMatrix<Number> & matrix,
1231 : DenseVector<Number> & rhs,
1232 : std::vector<dof_id_type> & elem_dofs,
1233 : bool asymmetric_constraint_rows = true) const;
1234 :
1235 : /**
1236 : * Constrains the element matrix and vector. This method requires
1237 : * the element matrix to be square, in which case the elem_dofs
1238 : * correspond to the global DOF indices of both the rows and
1239 : * columns of the element matrix. For this case the rows
1240 : * and columns of the matrix necessarily correspond to variables
1241 : * of the same approximation order.
1242 : *
1243 : * The heterogeneous version of this method creates linear systems in
1244 : * which heterogeneously constrained degrees of freedom will solve to
1245 : * their correct offset values, as would be appropriate for finding
1246 : * a solution to a linear problem in a single algebraic solve. The
1247 : * non-heterogeneous version of this method creates linear systems in
1248 : * which even heterogeneously constrained degrees of freedom are
1249 : * solved without offset values taken into account, as would be
1250 : * appropriate for finding linearized updates to a solution in which
1251 : * heterogeneous constraints are already satisfied.
1252 : *
1253 : * By default, the constraints for the primal solution of this
1254 : * system are used. If a non-negative \p qoi_index is passed in,
1255 : * then the constraints for the corresponding adjoint solution are
1256 : * used instead.
1257 : */
1258 : void heterogeneously_constrain_element_matrix_and_vector (DenseMatrix<Number> & matrix,
1259 : DenseVector<Number> & rhs,
1260 : std::vector<dof_id_type> & elem_dofs,
1261 : bool asymmetric_constraint_rows = true,
1262 : int qoi_index = -1) const;
1263 :
1264 : /*
1265 : * Backwards compatibility with misspelling.
1266 : */
1267 0 : void heterogenously_constrain_element_matrix_and_vector (DenseMatrix<Number> & matrix,
1268 : DenseVector<Number> & rhs,
1269 : std::vector<dof_id_type> & elem_dofs,
1270 : bool asymmetric_constraint_rows = true,
1271 : int qoi_index = -1) const
1272 : {
1273 : return this->heterogeneously_constrain_element_matrix_and_vector
1274 15527 : (matrix, rhs, elem_dofs, asymmetric_constraint_rows, qoi_index);
1275 : }
1276 :
1277 : /**
1278 : * Constrains the element vector. This method requires
1279 : * the element matrix to be square and not-yet-constrained, in which
1280 : * case the elem_dofs correspond to the global DOF indices of both
1281 : * the rows and columns of the element matrix.
1282 : *
1283 : * The heterogeneous version of this method creates linear systems in
1284 : * which heterogeneously constrained degrees of freedom will solve to
1285 : * their correct offset values, as would be appropriate for finding
1286 : * a solution to a linear problem in a single algebraic solve. The
1287 : * non-heterogeneous version of this method creates linear systems in
1288 : * which even heterogeneously constrained degrees of freedom are
1289 : * solved without offset values taken into account, as would be
1290 : * appropriate for finding linearized updates to a solution in which
1291 : * heterogeneous constraints are already satisfied.
1292 : *
1293 : * Note the sign difference from the nonlinear heterogeneous constraint
1294 : * method: Solving u:=K\f has the opposite sign convention from
1295 : * u:=u_in-J\r, and we apply heterogeneous constraints accordingly.
1296 : *
1297 : * By default, the constraints for the primal solution of this
1298 : * system are used. If a non-negative \p qoi_index is passed in,
1299 : * then the constraints for the corresponding adjoint solution are
1300 : * used instead.
1301 : */
1302 : void heterogeneously_constrain_element_vector (const DenseMatrix<Number> & matrix,
1303 : DenseVector<Number> & rhs,
1304 : std::vector<dof_id_type> & elem_dofs,
1305 : bool asymmetric_constraint_rows = true,
1306 : int qoi_index = -1) const;
1307 :
1308 : /*
1309 : * Backwards compatibility with misspelling.
1310 : */
1311 140 : void heterogenously_constrain_element_vector (const DenseMatrix<Number> & matrix,
1312 : DenseVector<Number> & rhs,
1313 : std::vector<dof_id_type> & elem_dofs,
1314 : bool asymmetric_constraint_rows = true,
1315 : int qoi_index = -1) const
1316 : {
1317 : return this->heterogeneously_constrain_element_vector
1318 1680 : (matrix, rhs, elem_dofs, asymmetric_constraint_rows, qoi_index);
1319 : }
1320 :
1321 : /**
1322 : * Constrains the element Jacobian and residual. The element
1323 : * Jacobian is square, and the elem_dofs should correspond to the
1324 : * global DOF indices of both the rows and columns of the element
1325 : * matrix.
1326 : *
1327 : * The residual-constraining version of this method creates linear
1328 : * systems in which heterogeneously constrained degrees of freedom
1329 : * create non-zero residual terms when not at their correct offset
1330 : * values, as would be appropriate for finding a solution to a
1331 : * nonlinear problem in a quasi-Newton solve.
1332 : *
1333 : * Note the sign difference from the linear heterogeneous constraint
1334 : * method: Solving u:=u_in-J\r has the opposite sign convention from
1335 : * u:=K\f, and we apply heterogeneous constraints accordingly.
1336 : *
1337 : * The \p solution vector passed in should be a serialized or
1338 : * ghosted primal solution
1339 : */
1340 : void heterogeneously_constrain_element_jacobian_and_residual (DenseMatrix<Number> & matrix,
1341 : DenseVector<Number> & rhs,
1342 : std::vector<dof_id_type> & elem_dofs,
1343 : NumericVector<Number> & solution_local) const;
1344 :
1345 : /**
1346 : * Constrains the element residual. The element Jacobian is square,
1347 : * and the elem_dofs should correspond to the global DOF indices of
1348 : * both the rows and columns of the element matrix.
1349 : *
1350 : * The residual-constraining version of this method creates linear
1351 : * systems in which heterogeneously constrained degrees of freedom
1352 : * create non-zero residual terms when not at their correct offset
1353 : * values, as would be appropriate for finding a solution to a
1354 : * nonlinear problem in a quasi-Newton solve.
1355 : *
1356 : * The \p solution vector passed in should be a serialized or
1357 : * ghosted primal solution
1358 : */
1359 : void heterogeneously_constrain_element_residual (DenseVector<Number> & rhs,
1360 : std::vector<dof_id_type> & elem_dofs,
1361 : NumericVector<Number> & solution_local) const;
1362 :
1363 :
1364 : /**
1365 : * Constrains the element residual. The element Jacobian is square,
1366 : * and the elem_dofs should correspond to the global DOF indices of
1367 : * both the rows and columns of the element matrix, and the dof
1368 : * constraint should not include any heterogeneous terms.
1369 : *
1370 : * The residual-constraining version of this method creates linear
1371 : * systems in which heterogeneously constrained degrees of freedom
1372 : * create non-zero residual terms when not at their correct offset
1373 : * values, as would be appropriate for finding a solution to a
1374 : * nonlinear problem in a quasi-Newton solve.
1375 : *
1376 : * The \p solution vector passed in should be a serialized or
1377 : * ghosted primal solution
1378 : */
1379 : void constrain_element_residual (DenseVector<Number> & rhs,
1380 : std::vector<dof_id_type> & elem_dofs,
1381 : NumericVector<Number> & solution_local) const;
1382 :
1383 : /**
1384 : * Constrains a dyadic element matrix B = v w'. This method
1385 : * requires the element matrix to be square, in which case the
1386 : * elem_dofs correspond to the global DOF indices of both the rows
1387 : * and columns of the element matrix. For this case the rows and
1388 : * columns of the matrix necessarily correspond to variables of the
1389 : * same approximation order.
1390 : */
1391 : void constrain_element_dyad_matrix (DenseVector<Number> & v,
1392 : DenseVector<Number> & w,
1393 : std::vector<dof_id_type> & row_dofs,
1394 : bool asymmetric_constraint_rows = true) const;
1395 :
1396 : /**
1397 : * Does not actually constrain anything, but modifies \p dofs in the
1398 : * same way as any of the constrain functions would do, i.e. adds
1399 : * those dofs in terms of which any of the existing dofs is
1400 : * constrained.
1401 : */
1402 : void constrain_nothing (std::vector<dof_id_type> & dofs) const;
1403 :
1404 : /**
1405 : * Constrains the numeric vector \p v, which represents a solution defined on
1406 : * the mesh. This may need to be used after a linear solve, if your linear
1407 : * solver's solutions do not satisfy your DoF constraints to a tight enough
1408 : * tolerance.
1409 : *
1410 : * If \p v == nullptr, the system solution vector is constrained
1411 : *
1412 : * If \p homogeneous == true, heterogeneous constraints are enforced
1413 : * as if they were homogeneous. This might be appropriate for e.g. a
1414 : * vector representing a difference between two
1415 : * heterogeneously-constrained solutions.
1416 : */
1417 : void enforce_constraints_exactly (const System & system,
1418 : NumericVector<Number> * v = nullptr,
1419 : bool homogeneous = false) const;
1420 :
1421 : /**
1422 : * Heterogeneously constrains the numeric vector \p v, which
1423 : * represents an adjoint solution defined on the mesh for quantity
1424 : * fo interest \p q. For homogeneous constraints, use \p
1425 : * enforce_constraints_exactly instead
1426 : */
1427 : void enforce_adjoint_constraints_exactly (NumericVector<Number> & v,
1428 : unsigned int q) const;
1429 :
1430 : void enforce_constraints_on_residual (const NonlinearImplicitSystem & system,
1431 : NumericVector<Number> * rhs,
1432 : NumericVector<Number> const * solution,
1433 : bool homogeneous = true) const;
1434 :
1435 : void enforce_constraints_on_jacobian (const NonlinearImplicitSystem & system,
1436 : SparseMatrix<Number> * jac) const;
1437 :
1438 : #ifdef LIBMESH_ENABLE_PERIODIC
1439 :
1440 : //--------------------------------------------------------------------
1441 : // PeriodicBoundary-specific methods
1442 :
1443 : /**
1444 : * Adds a copy of the specified periodic boundary to the system.
1445 : */
1446 : void add_periodic_boundary (const PeriodicBoundaryBase & periodic_boundary);
1447 :
1448 : /**
1449 : * Add a periodic boundary pair
1450 : *
1451 : * \param boundary - primary boundary
1452 : * \param inverse_boundary - inverse boundary
1453 : */
1454 : void add_periodic_boundary (const PeriodicBoundaryBase & boundary, const PeriodicBoundaryBase & inverse_boundary);
1455 :
1456 : /**
1457 : * \returns \p true if the boundary given by \p boundaryid is periodic,
1458 : * false otherwise
1459 : */
1460 : bool is_periodic_boundary (const boundary_id_type boundaryid) const;
1461 :
1462 : PeriodicBoundaries * get_periodic_boundaries()
1463 : {
1464 : return _periodic_boundaries.get();
1465 : }
1466 :
1467 : const PeriodicBoundaries * get_periodic_boundaries() const
1468 : {
1469 : return _periodic_boundaries.get();
1470 : }
1471 :
1472 : #endif // LIBMESH_ENABLE_PERIODIC
1473 :
1474 :
1475 : #ifdef LIBMESH_ENABLE_DIRICHLET
1476 :
1477 : //--------------------------------------------------------------------
1478 : // DirichletBoundary-specific methods
1479 :
1480 : /**
1481 : * Adds a copy of the specified Dirichlet boundary to the system.
1482 : *
1483 : * The constraints implied by DirichletBoundary objects are imposed
1484 : * in the same order in which DirichletBoundary objects are added to
1485 : * the DofMap. When multiple DirichletBoundary objects would impose
1486 : * competing constraints on a given DOF, the *first*
1487 : * DirichletBoundary to constrain the DOF "wins". This distinction
1488 : * is important when e.g. two surfaces (sidesets) intersect. The
1489 : * nodes on the intersection will be constrained according to
1490 : * whichever sideset's DirichletBoundary object was added to the
1491 : * DofMap first.
1492 : */
1493 : void add_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary);
1494 :
1495 : /**
1496 : * Adds a copy of the specified Dirichlet boundary to the system,
1497 : * corresponding to the adjoint problem defined by Quantity of
1498 : * Interest \p q.
1499 : */
1500 : void add_adjoint_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary,
1501 : unsigned int q);
1502 :
1503 : /**
1504 : * Removes the specified Dirichlet boundary from the system.
1505 : */
1506 : void remove_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary);
1507 :
1508 : /**
1509 : * Removes from the system the specified Dirichlet boundary for the
1510 : * adjoint equation defined by Quantity of interest index q
1511 : */
1512 : void remove_adjoint_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary,
1513 : unsigned int q);
1514 :
1515 : const DirichletBoundaries * get_dirichlet_boundaries() const
1516 : {
1517 : return _dirichlet_boundaries.get();
1518 : }
1519 :
1520 24 : DirichletBoundaries * get_dirichlet_boundaries()
1521 : {
1522 24 : return _dirichlet_boundaries.get();
1523 : }
1524 :
1525 : bool has_adjoint_dirichlet_boundaries(unsigned int q) const;
1526 :
1527 : const DirichletBoundaries *
1528 : get_adjoint_dirichlet_boundaries(unsigned int q) const;
1529 :
1530 : DirichletBoundaries *
1531 : get_adjoint_dirichlet_boundaries(unsigned int q);
1532 :
1533 : /**
1534 : * Check that all the ids in dirichlet_bcids are actually present in the mesh.
1535 : * If not, this will throw an error.
1536 : */
1537 : void check_dirichlet_bcid_consistency (const MeshBase & mesh,
1538 : const DirichletBoundary & boundary) const;
1539 : #endif // LIBMESH_ENABLE_DIRICHLET
1540 :
1541 :
1542 : #ifdef LIBMESH_ENABLE_AMR
1543 :
1544 : //--------------------------------------------------------------------
1545 : // AMR-specific methods
1546 :
1547 : /**
1548 : * After a mesh is refined and repartitioned it is possible that the
1549 : * \p _send_list will need to be augmented. This is the case when an
1550 : * element is refined and its children end up on different processors
1551 : * than the parent. These children will need values from the parent
1552 : * when projecting the solution onto the refined mesh, hence the parent's
1553 : * DOF indices need to be included in the \p _send_list.
1554 : */
1555 : // void augment_send_list_for_projection(const MeshBase &);
1556 :
1557 : #ifdef LIBMESH_ENABLE_AMR
1558 :
1559 : /**
1560 : * Fills the vector di with the global degree of freedom indices
1561 : * for the element using the \p DofMap::old_dof_object.
1562 : * If no variable number is specified then all
1563 : * variables are returned.
1564 : */
1565 : void old_dof_indices (const Elem * const elem,
1566 : std::vector<dof_id_type> & di,
1567 : const unsigned int vn = libMesh::invalid_uint) const;
1568 :
1569 : #endif // LIBMESH_ENABLE_AMR
1570 :
1571 : /**
1572 : * Constrains degrees of freedom on side \p s of element \p elem which
1573 : * correspond to variable number \p var and to p refinement levels
1574 : * above \p p.
1575 : */
1576 : void constrain_p_dofs (unsigned int var,
1577 : const Elem * elem,
1578 : unsigned int s,
1579 : unsigned int p);
1580 :
1581 : #endif // LIBMESH_ENABLE_AMR
1582 :
1583 : /**
1584 : * Reinitialize the underlying data structures conformal to the current mesh.
1585 : */
1586 : void reinit
1587 : (MeshBase & mesh,
1588 : const std::map<const Node *, std::set<subdomain_id_type>> &
1589 : constraining_subdomains);
1590 :
1591 : /**
1592 : * Free all new memory associated with the object, but restore its
1593 : * original state, with the mesh pointer and any default ghosting.
1594 : */
1595 : virtual void clear () override;
1596 :
1597 : /**
1598 : * Prints summary info about the sparsity bandwidth and constraints.
1599 : */
1600 : void print_info(std::ostream & os=libMesh::out) const;
1601 :
1602 : /**
1603 : * Gets summary info about the sparsity bandwidth and constraints.
1604 : */
1605 : std::string get_info() const;
1606 :
1607 : /**
1608 : * Degree of freedom coupling. If left empty each DOF
1609 : * couples to all others. Can be used to reduce memory
1610 : * requirements for sparse matrices. DOF 0 might only
1611 : * couple to itself, in which case \p dof_coupling(0,0)
1612 : * should be 1 and \p dof_coupling(0,j) = 0 for j not equal
1613 : * to 0.
1614 : *
1615 : * This variable is named as though it were class private,
1616 : * but it is in the public interface. Also there are no
1617 : * public methods for accessing it... This typically means
1618 : * you should only use it if you know what you are doing.
1619 : */
1620 : CouplingMatrix * _dof_coupling;
1621 :
1622 : /**
1623 : * \returns The number of the system we are responsible for.
1624 : */
1625 : unsigned int sys_number() const;
1626 :
1627 : /**
1628 : * Builds a sparsity pattern for matrices using the current
1629 : * degree-of-freedom numbering and coupling.
1630 : *
1631 : * By default, ignores constraint equations, for build speed; this
1632 : * is valid for the combination of !need_full_sparsity_pattern and
1633 : * constraints which only come from periodic boundary conditions and
1634 : * adaptive mesh refinement, where matrix constraint adds some
1635 : * matrix entries but removes equally many (or more) other entries.
1636 : *
1637 : * Can be told to calculate sparsity for the constrained matrix,
1638 : * which may be necessary in the case of spline control node
1639 : * constraints or sufficiently many user constraints.
1640 : */
1641 : std::unique_ptr<SparsityPattern::Build> build_sparsity(const MeshBase & mesh,
1642 : bool calculate_constrained = false,
1643 : bool use_condensed_system = false) const;
1644 :
1645 : /**
1646 : * Describe whether the given variable group should be p-refined. If this API is not called with
1647 : * \p false, the default is to p-refine
1648 : */
1649 : void should_p_refine(unsigned int g, bool p_refine);
1650 :
1651 : /**
1652 : * Whether the given variable group should be p-refined
1653 : */
1654 : bool should_p_refine(unsigned int g) const;
1655 :
1656 : /**
1657 : * Whether the given variable should be p-refined
1658 : */
1659 : bool should_p_refine_var(unsigned int var) const;
1660 :
1661 : // Prevent bad user implicit conversions
1662 : void should_p_refine(FEFamily, bool) = delete;
1663 : void should_p_refine(Order, bool) = delete;
1664 : bool should_p_refine(FEFamily) const = delete;
1665 : bool should_p_refine(Order) const = delete;
1666 :
1667 : /**
1668 : * Add a static condensation class
1669 : */
1670 : void create_static_condensation(MeshBase & mesh, System & system);
1671 :
1672 : /**
1673 : * Checks whether we have static condensation
1674 : */
1675 50582 : bool has_static_condensation() const { return _sc.get(); }
1676 :
1677 : /**
1678 : * @returns the static condensation class. This should have been already added with a call to \p
1679 : * add_static_condensation()
1680 : */
1681 : StaticCondensationDofMap & get_static_condensation();
1682 :
1683 : /**
1684 : * Calls reinit on the static condensation map if it exists
1685 : */
1686 : void reinit_static_condensation();
1687 :
1688 : private:
1689 :
1690 : /**
1691 : * Helper function that gets the dof indices on the current element
1692 : * for a non-SCALAR type variable, where the variable is identified
1693 : * by its variable group number \p vg and its offset \p vig from the
1694 : * first variable in that group.
1695 : *
1696 : * In DEBUG mode, the tot_size parameter will add up the total
1697 : * number of dof indices that should have been added to di, and v
1698 : * will be the variable number corresponding to vg and vig.
1699 : */
1700 : void _dof_indices (const Elem & elem,
1701 : int p_level,
1702 : std::vector<dof_id_type> & di,
1703 : const unsigned int vg,
1704 : const unsigned int vig,
1705 : const Node * const * nodes,
1706 : unsigned int n_nodes,
1707 : const unsigned int v
1708 : #ifdef DEBUG
1709 : ,
1710 : std::size_t & tot_size
1711 : #endif
1712 : ) const;
1713 :
1714 : /**
1715 : * As above except a \p field_dofs_functor must be provided. This method is useful when the caller
1716 : * wants to do more than simply fill a degree of freedom container
1717 : * @param field_dofs_functor This functor has the interface:
1718 : * void field_dofs_functor(const Elem & elem,
1719 : * const unsigned int node_num,
1720 : * const unsigned int var_num,
1721 : * std::vector<dof_id_type> & di,
1722 : * const dof_id_type field_dof)
1723 : * where \p field_dof represents a field degree of freedom to act on and
1724 : * is associated with \p node_num and \p var_num. If the degree of
1725 : * freedom is elemental than \p node_num will be \p invalid_uint. \p di
1726 : * is the degree of freedom container provided to the \p _dof_indices
1727 : * method
1728 : */
1729 : template <typename FieldDofsFunctor>
1730 : void _dof_indices (const Elem & elem,
1731 : int p_level,
1732 : std::vector<dof_id_type> & di,
1733 : const unsigned int vg,
1734 : const unsigned int vig,
1735 : const Node * const * nodes,
1736 : unsigned int n_nodes,
1737 : const unsigned int v,
1738 : #ifdef DEBUG
1739 : std::size_t & tot_size,
1740 : #endif
1741 : FieldDofsFunctor field_dofs_functor) const;
1742 :
1743 : /**
1744 : * Helper function that implements the element-nodal versions of
1745 : * dof_indices and old_dof_indices
1746 : */
1747 : void _node_dof_indices (const Elem & elem,
1748 : unsigned int n,
1749 : const DofObject & obj,
1750 : std::vector<dof_id_type> & di,
1751 : const unsigned int vn) const;
1752 :
1753 : /**
1754 : * Invalidates all active DofObject dofs for this system
1755 : */
1756 : void invalidate_dofs(MeshBase & mesh) const;
1757 :
1758 : /**
1759 : * \returns The Node pointer with index \p i from the \p mesh.
1760 : */
1761 : DofObject * node_ptr(MeshBase & mesh, dof_id_type i) const;
1762 :
1763 : /**
1764 : * \returns The Elem pointer with index \p i from the \p mesh.
1765 : */
1766 : DofObject * elem_ptr(MeshBase & mesh, dof_id_type i) const;
1767 :
1768 : /**
1769 : * A member function type like \p node_ptr() or \p elem_ptr().
1770 : */
1771 : typedef DofObject * (DofMap::*dofobject_accessor)
1772 : (MeshBase & mesh, dof_id_type i) const;
1773 :
1774 : /**
1775 : * Helper function for distributing dofs in parallel
1776 : */
1777 : template<typename iterator_type>
1778 : void set_nonlocal_dof_objects(iterator_type objects_begin,
1779 : iterator_type objects_end,
1780 : MeshBase & mesh,
1781 : dofobject_accessor objects);
1782 :
1783 : /**
1784 : * We may have mesh constraint rows with dependent nodes in one
1785 : * subdomain but dependency nodes in another subdomain, and we may
1786 : * have variables whose subdomain restriction includes the dependent
1787 : * subdomain but not the dependency. In those cases we need to
1788 : * place degrees of freedom on dependency nodes anyway.
1789 : *
1790 : * The set value for node n will include all subdomain ids of
1791 : * elements with nodes in subdomains constrained by n.
1792 : *
1793 : * We use a map<set> rather than a multimap here because we expect
1794 : * to be inserting the same subdomain multiple times and we don't
1795 : * need duplicate values.
1796 : */
1797 : std::map<const Node *, std::set<subdomain_id_type>>
1798 : calculate_constraining_subdomains();
1799 :
1800 : /**
1801 : * Distributes the global degrees of freedom, for dofs on
1802 : * this processor. In this format the local
1803 : * degrees of freedom are in a contiguous block for each
1804 : * variable in the system.
1805 : * Starts at index next_free_dof, and increments it to
1806 : * the post-final index.
1807 : *
1808 : * Uses the provided constraining_subdomains map from
1809 : * calculate_constraining_subdomains() to ensure allocation of all
1810 : * DoFs on constraining nodes.
1811 : */
1812 : void distribute_local_dofs_var_major
1813 : (dof_id_type & next_free_dof,
1814 : MeshBase & mesh,
1815 : const std::map<const Node *, std::set<subdomain_id_type>> &
1816 : constraining_subdomains);
1817 :
1818 : /**
1819 : * Distributes the global degrees of freedom for dofs on this
1820 : * processor. In this format all the degrees of freedom at a
1821 : * node/element are in contiguous blocks. Starts at index \p
1822 : * next_free_dof, and increments it to the post-final index. If \p
1823 : * build_send_list is \p true, builds the send list. If \p false,
1824 : * clears and reserves the send list.
1825 : *
1826 : * Uses the provided constraining_subdomains map from
1827 : * calculate_constraining_subdomains() to ensure allocation of all
1828 : * DoFs on constraining nodes.
1829 : *
1830 : * \note The degrees of freedom for a given variable are not in
1831 : * contiguous blocks, as in the case of \p distribute_local_dofs_var_major.
1832 : */
1833 : void distribute_local_dofs_node_major
1834 : (dof_id_type & next_free_dof,
1835 : MeshBase & mesh,
1836 : const std::map<const Node *, std::set<subdomain_id_type>> &
1837 : constraining_subdomains);
1838 :
1839 : /*
1840 : * Helper method for the above two to count + distriubte SCALAR dofs
1841 : */
1842 : void distribute_scalar_dofs (dof_id_type & next_free_dof);
1843 :
1844 : #ifdef DEBUG
1845 : /*
1846 : * Internal assertions for distribute_local_dofs_*
1847 : */
1848 : void assert_no_nodes_missed(MeshBase & mesh);
1849 : #endif
1850 :
1851 : /*
1852 : * A utility method for obtaining a set of elements to ghost along
1853 : * with merged coupling matrices.
1854 : */
1855 : typedef std::set<std::unique_ptr<CouplingMatrix>, Utility::CompareUnderlying> CouplingMatricesSet;
1856 : static void
1857 : merge_ghost_functor_outputs (GhostingFunctor::map_type & elements_to_ghost,
1858 : CouplingMatricesSet & temporary_coupling_matrices,
1859 : const GhostingFunctorIterator & gf_begin,
1860 : const GhostingFunctorIterator & gf_end,
1861 : const MeshBase::const_element_iterator & elems_begin,
1862 : const MeshBase::const_element_iterator & elems_end,
1863 : processor_id_type p);
1864 :
1865 : /**
1866 : * Adds entries to the \p _send_list vector corresponding to DoFs
1867 : * on elements neighboring the current processor.
1868 : */
1869 : void add_neighbors_to_send_list(MeshBase & mesh);
1870 :
1871 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
1872 :
1873 : /**
1874 : * Build the constraint matrix C associated with the element
1875 : * degree of freedom indices elem_dofs. The optional parameter
1876 : * \p called_recursively should be left at the default value
1877 : * \p false. This is used to handle the special case of
1878 : * an element's degrees of freedom being constrained in terms
1879 : * of other, local degrees of freedom. The usual case is
1880 : * for an elements DOFs to be constrained by some other,
1881 : * external DOFs.
1882 : */
1883 : void build_constraint_matrix (DenseMatrix<Number> & C,
1884 : std::vector<dof_id_type> & elem_dofs,
1885 : const bool called_recursively=false) const;
1886 :
1887 : /**
1888 : * Build the constraint matrix C and the forcing vector H
1889 : * associated with the element degree of freedom indices elem_dofs.
1890 : * The optional parameter \p called_recursively should be left at
1891 : * the default value \p false. This is used to handle the special
1892 : * case of an element's degrees of freedom being constrained in
1893 : * terms of other, local degrees of freedom. The usual case is for
1894 : * an elements DOFs to be constrained by some other, external DOFs
1895 : * and/or Dirichlet conditions.
1896 : *
1897 : * The forcing vector will depend on which solution's heterogeneous
1898 : * constraints are being applied. For the default \p qoi_index this
1899 : * will be the primal solution; for \p qoi_index >= 0 the
1900 : * corresponding adjoint solution's constraints will be used.
1901 : */
1902 : void build_constraint_matrix_and_vector (DenseMatrix<Number> & C,
1903 : DenseVector<Number> & H,
1904 : std::vector<dof_id_type> & elem_dofs,
1905 : int qoi_index = -1,
1906 : const bool called_recursively=false) const;
1907 :
1908 : /**
1909 : * Finds all the DOFS associated with the element DOFs elem_dofs.
1910 : * This will account for off-element couplings via hanging nodes.
1911 : */
1912 : void find_connected_dofs (std::vector<dof_id_type> & elem_dofs) const;
1913 :
1914 : /**
1915 : * Finds all the DofObjects associated with the set in \p objs.
1916 : * This will account for off-element couplings via hanging nodes.
1917 : */
1918 : void find_connected_dof_objects (std::vector<const DofObject *> & objs) const;
1919 :
1920 : /**
1921 : * Adds entries to the \p _send_list vector corresponding to DoFs
1922 : * which are dependencies for constraint equations on the current
1923 : * processor.
1924 : */
1925 : void add_constraints_to_send_list();
1926 :
1927 : /**
1928 : * Adds any spline constraints from the Mesh to our DoF constraints.
1929 : * If any Dirichlet constraints exist on spline-constrained nodes,
1930 : * l2-projects those constraints onto the spline basis.
1931 : */
1932 : void process_mesh_constraint_rows(const MeshBase & mesh);
1933 :
1934 : #endif // LIBMESH_ENABLE_CONSTRAINTS
1935 :
1936 : /**
1937 : * This flag indicates whether or not we do an opt-mode check for
1938 : * the presence of constraint loops, i.e. cases where the constraint
1939 : * graph is cyclic.
1940 : */
1941 : bool _error_on_constraint_loop;
1942 :
1943 : /**
1944 : * This flag indicates whether or not we explicitly take constraint
1945 : * equations into account when computing a sparsity pattern.
1946 : */
1947 : bool _constrained_sparsity_construction;
1948 :
1949 : /**
1950 : * The finite element type for each variable.
1951 : */
1952 : std::vector<Variable> _variables;
1953 :
1954 : /**
1955 : * The finite element type for each variable group.
1956 : */
1957 : std::vector<VariableGroup> _variable_groups;
1958 :
1959 : /**
1960 : * The variable group number for each variable.
1961 : */
1962 : std::vector<unsigned int> _variable_group_numbers;
1963 :
1964 : /**
1965 : * A map from variable number to variable group number
1966 : */
1967 : std::unordered_map<unsigned int, unsigned int> _var_to_vg;
1968 :
1969 : /**
1970 : * The number of the system we manage DOFs for.
1971 : */
1972 : const unsigned int _sys_number;
1973 :
1974 : /**
1975 : * The mesh that system uses.
1976 : */
1977 : MeshBase & _mesh;
1978 :
1979 : /**
1980 : * Additional matrices handled by this object. These pointers do \e
1981 : * not handle the memory, instead, \p System, who
1982 : * told \p DofMap about them, owns them.
1983 : */
1984 : std::vector<SparseMatrix<Number> * > _matrices;
1985 :
1986 : /**
1987 : * First DOF index for SCALAR variable v, or garbage for non-SCALAR
1988 : * variable v
1989 : */
1990 : std::vector<dof_id_type> _first_scalar_df;
1991 :
1992 : /**
1993 : * A list containing all the global DOF indices that affect the
1994 : * solution on my processor.
1995 : */
1996 : std::vector<dof_id_type> _send_list;
1997 :
1998 : /**
1999 : * Function object to call to add extra entries to the sparsity pattern
2000 : */
2001 : SparsityPattern::AugmentSparsityPattern * _augment_sparsity_pattern;
2002 :
2003 : /**
2004 : * A function pointer to a function to call to add extra entries to the sparsity pattern
2005 : */
2006 : void (*_extra_sparsity_function)(SparsityPattern::Graph &,
2007 : std::vector<dof_id_type> & n_nz,
2008 : std::vector<dof_id_type> & n_oz,
2009 : void *);
2010 : /**
2011 : * A pointer associated with the extra sparsity that can optionally be passed in
2012 : */
2013 : void * _extra_sparsity_context;
2014 :
2015 : /**
2016 : * Function object to call to add extra entries to the send list
2017 : */
2018 : AugmentSendList * _augment_send_list;
2019 :
2020 : /**
2021 : * A function pointer to a function to call to add extra entries to the send list
2022 : */
2023 : void (*_extra_send_list_function)(std::vector<dof_id_type> &, void *);
2024 :
2025 : /**
2026 : * A pointer associated with the extra send list that can optionally be passed in
2027 : */
2028 : void * _extra_send_list_context;
2029 :
2030 : /**
2031 : * The default coupling GhostingFunctor, used to implement standard
2032 : * libMesh sparsity pattern construction.
2033 : *
2034 : * We use a std::unique_ptr here to reduce header dependencies.
2035 : */
2036 : std::unique_ptr<DefaultCoupling> _default_coupling;
2037 :
2038 : /**
2039 : * The default algebraic GhostingFunctor, used to implement standard
2040 : * libMesh send_list construction.
2041 : *
2042 : * We use a std::unique_ptr here to reduce header dependencies.
2043 : */
2044 : std::unique_ptr<DefaultCoupling> _default_evaluating;
2045 :
2046 : /**
2047 : * The list of all GhostingFunctor objects to be used when
2048 : * distributing ghosted vectors.
2049 : *
2050 : * The library should automatically refer these functors to the
2051 : * MeshBase, too, so any algebraically ghosted dofs will live on
2052 : * geometrically ghosted elements.
2053 : *
2054 : * Keep these in a vector so any parallel computation is done in the
2055 : * same order on all processors.
2056 : */
2057 : std::vector<GhostingFunctor *> _algebraic_ghosting_functors;
2058 :
2059 : /**
2060 : * The list of all GhostingFunctor objects to be used when
2061 : * coupling degrees of freedom in matrix sparsity patterns.
2062 : *
2063 : * These objects will *also* be used as algebraic ghosting functors,
2064 : * but not vice-versa.
2065 : *
2066 : * The library should automatically refer these functors to the
2067 : * MeshBase, too, so any dofs coupled to local dofs will live on
2068 : * geometrically ghosted elements.
2069 : */
2070 : std::vector<GhostingFunctor *> _coupling_functors;
2071 :
2072 : /**
2073 : * Hang on to references to any GhostingFunctor objects we were
2074 : * passed in shared_ptr form
2075 : */
2076 : std::map<GhostingFunctor *, std::shared_ptr<GhostingFunctor> > _shared_functors;
2077 :
2078 : /**
2079 : * Default false; set to true if any attached matrix requires a full
2080 : * sparsity pattern.
2081 : */
2082 : bool need_full_sparsity_pattern;
2083 :
2084 : /**
2085 : * The sparsity pattern of the global matrix. If
2086 : * need_full_sparsity_pattern is true, we save the entire sparse
2087 : * graph here. Otherwise we save just the n_nz and n_oz vectors.
2088 : */
2089 : std::unique_ptr<SparsityPattern::Build> _sp;
2090 :
2091 : /**
2092 : * The total number of SCALAR dofs associated to
2093 : * all SCALAR variables.
2094 : */
2095 : dof_id_type _n_SCALAR_dofs;
2096 :
2097 : #ifdef LIBMESH_ENABLE_AMR
2098 :
2099 : /**
2100 : * First old DOF index for SCALAR variable v, or garbage for
2101 : * non-SCALAR variable v
2102 : */
2103 : std::vector<dof_id_type> _first_old_scalar_df;
2104 :
2105 : /**
2106 : * A container of variable groups that we should not p-refine
2107 : */
2108 : std::unordered_set<unsigned int> _dont_p_refine;
2109 : #endif
2110 :
2111 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
2112 : /**
2113 : * Data structure containing DOF constraints. The ith
2114 : * entry is the constraint matrix row for DOF i.
2115 : */
2116 : DofConstraints _dof_constraints, _stashed_dof_constraints;
2117 :
2118 : DofConstraintValueMap _primal_constraint_values;
2119 :
2120 : AdjointDofConstraintValues _adjoint_constraint_values;
2121 : #endif
2122 :
2123 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2124 : /**
2125 : * Data structure containing DofObject constraints.
2126 : */
2127 : NodeConstraints _node_constraints;
2128 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2129 :
2130 :
2131 : #ifdef LIBMESH_ENABLE_PERIODIC
2132 : /**
2133 : * Data structure containing periodic boundaries. The ith
2134 : * entry is the constraint matrix row for boundaryid i.
2135 : */
2136 : std::unique_ptr<PeriodicBoundaries> _periodic_boundaries;
2137 : #endif
2138 :
2139 : #ifdef LIBMESH_ENABLE_DIRICHLET
2140 : /**
2141 : * Data structure containing Dirichlet functions. The ith
2142 : * entry is the constraint matrix row for boundaryid i.
2143 : */
2144 : std::unique_ptr<DirichletBoundaries> _dirichlet_boundaries;
2145 :
2146 : /**
2147 : * Data structure containing Dirichlet functions. The ith
2148 : * entry is the constraint matrix row for boundaryid i.
2149 : */
2150 : std::vector<std::unique_ptr<DirichletBoundaries>> _adjoint_dirichlet_boundaries;
2151 : #endif
2152 :
2153 : friend class SparsityPattern::Build;
2154 :
2155 : /**
2156 : * Bools to indicate if we override the --implicit_neighbor_dofs
2157 : * commandline options.
2158 : */
2159 : bool _implicit_neighbor_dofs_initialized;
2160 : bool _implicit_neighbor_dofs;
2161 :
2162 : /**
2163 : * Flag which determines whether we should do some additional
2164 : * checking of the consistency of the DirichletBoundary objects
2165 : * added by the user. Defaults to true, but can be disabled in cases
2166 : * where you only want to add DirichletBoundary objects "locally"
2167 : * and can guarantee that no repartitioning will be done, since
2168 : * repartitioning could cause processors to own new boundary sides
2169 : * for which they no longer have the proper DirichletBoundary
2170 : * objects stored.
2171 : */
2172 : bool _verify_dirichlet_bc_consistency;
2173 :
2174 : /// Static condensation class
2175 : std::unique_ptr<StaticCondensationDofMap> _sc;
2176 : };
2177 :
2178 :
2179 : // ------------------------------------------------------------
2180 : // Dof Map inline member functions
2181 : inline
2182 42679894 : unsigned int DofMap::sys_number() const
2183 : {
2184 293042526 : return _sys_number;
2185 : }
2186 :
2187 :
2188 :
2189 : inline
2190 47455298 : const VariableGroup & DofMap::variable_group (const unsigned int g) const
2191 : {
2192 47455298 : libmesh_assert_less (g, _variable_groups.size());
2193 :
2194 302541438 : return _variable_groups[g];
2195 : }
2196 :
2197 :
2198 :
2199 : inline
2200 14669140 : const Variable & DofMap::variable (const unsigned int c) const
2201 : {
2202 1006624 : libmesh_assert_less (c, _variables.size());
2203 :
2204 15498471 : return _variables[c];
2205 : }
2206 :
2207 :
2208 :
2209 : inline
2210 : Order DofMap::variable_order (const unsigned int c) const
2211 : {
2212 : libmesh_assert_less (c, _variables.size());
2213 :
2214 : return _variables[c].type().order;
2215 : }
2216 :
2217 :
2218 :
2219 : inline
2220 : Order DofMap::variable_group_order (const unsigned int vg) const
2221 : {
2222 : libmesh_assert_less (vg, _variable_groups.size());
2223 :
2224 : return _variable_groups[vg].type().order;
2225 : }
2226 :
2227 :
2228 :
2229 : inline
2230 519326 : const FEType & DofMap::variable_type (const unsigned int c) const
2231 : {
2232 519326 : libmesh_assert_less (c, _variables.size());
2233 :
2234 27341733 : return _variables[c].type();
2235 : }
2236 :
2237 :
2238 :
2239 : inline
2240 : const FEType & DofMap::variable_group_type (const unsigned int vg) const
2241 : {
2242 : libmesh_assert_less (vg, _variable_groups.size());
2243 :
2244 : return _variable_groups[vg].type();
2245 : }
2246 :
2247 :
2248 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
2249 :
2250 :
2251 : inline
2252 1181091 : bool DofMap::is_constrained_node (const Node *
2253 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2254 : node
2255 : #endif
2256 : ) const
2257 : {
2258 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2259 1181091 : if (_node_constraints.count(node))
2260 13045 : return true;
2261 : #endif
2262 :
2263 1168046 : return false;
2264 : }
2265 :
2266 :
2267 : inline
2268 50803408 : bool DofMap::is_constrained_dof (const dof_id_type dof) const
2269 : {
2270 50803408 : if (_dof_constraints.count(dof))
2271 5418954 : return true;
2272 :
2273 45384454 : return false;
2274 : }
2275 :
2276 :
2277 : inline
2278 83110 : bool DofMap::has_heterogeneous_adjoint_constraints (const unsigned int qoi_num) const
2279 : {
2280 : AdjointDofConstraintValues::const_iterator it =
2281 4796 : _adjoint_constraint_values.find(qoi_num);
2282 87906 : if (it == _adjoint_constraint_values.end())
2283 1892 : return false;
2284 29012 : if (it->second.empty())
2285 25920 : return false;
2286 :
2287 24 : return true;
2288 : }
2289 :
2290 :
2291 : inline
2292 1415182 : Number DofMap::has_heterogeneous_adjoint_constraint (const unsigned int qoi_num,
2293 : const dof_id_type dof) const
2294 : {
2295 : AdjointDofConstraintValues::const_iterator it =
2296 1094202 : _adjoint_constraint_values.find(qoi_num);
2297 1415182 : if (it != _adjoint_constraint_values.end())
2298 : {
2299 : DofConstraintValueMap::const_iterator rhsit =
2300 104428 : it->second.find(dof);
2301 425408 : if (rhsit == it->second.end())
2302 154302 : return 0;
2303 : else
2304 1820 : return rhsit->second;
2305 : }
2306 :
2307 989774 : return 0;
2308 : }
2309 :
2310 :
2311 :
2312 : inline
2313 : DofConstraintValueMap & DofMap::get_primal_constraint_values()
2314 : {
2315 : return _primal_constraint_values;
2316 : }
2317 :
2318 :
2319 :
2320 : #else
2321 :
2322 : //--------------------------------------------------------------------
2323 : // Constraint-specific methods get inlined into nothing if
2324 : // constraints are disabled, so there's no reason for users not to
2325 : // use them.
2326 :
2327 : inline void DofMap::constrain_element_matrix (DenseMatrix<Number> &,
2328 : std::vector<dof_id_type> &,
2329 : bool) const {}
2330 :
2331 : inline void DofMap::constrain_element_matrix (DenseMatrix<Number> &,
2332 : std::vector<dof_id_type> &,
2333 : std::vector<dof_id_type> &,
2334 : bool) const {}
2335 :
2336 : inline void DofMap::constrain_element_vector (DenseVector<Number> &,
2337 : std::vector<dof_id_type> &,
2338 : bool) const {}
2339 :
2340 : inline void DofMap::constrain_element_matrix_and_vector (DenseMatrix<Number> &,
2341 : DenseVector<Number> &,
2342 : std::vector<dof_id_type> &,
2343 : bool) const {}
2344 :
2345 : inline void DofMap::heterogeneously_constrain_element_matrix_and_vector
2346 : (DenseMatrix<Number> &, DenseVector<Number> &,
2347 : std::vector<dof_id_type> &, bool, int) const {}
2348 :
2349 : inline void DofMap::heterogeneously_constrain_element_vector
2350 : (const DenseMatrix<Number> &, DenseVector<Number> &,
2351 : std::vector<dof_id_type> &, bool, int) const {}
2352 :
2353 : inline void DofMap::constrain_element_dyad_matrix (DenseVector<Number> &,
2354 : DenseVector<Number> &,
2355 : std::vector<dof_id_type> &,
2356 : bool) const {}
2357 :
2358 : inline void DofMap::constrain_nothing (std::vector<dof_id_type> &) const {}
2359 :
2360 : inline void DofMap::enforce_constraints_exactly (const System &,
2361 : NumericVector<Number> *,
2362 : bool) const {}
2363 :
2364 : inline void DofMap::enforce_adjoint_constraints_exactly (NumericVector<Number> &,
2365 : unsigned int) const {}
2366 :
2367 :
2368 : inline void DofMap::enforce_constraints_on_residual
2369 : (const NonlinearImplicitSystem &,
2370 : NumericVector<Number> *,
2371 : NumericVector<Number> const *,
2372 : bool) const {}
2373 :
2374 : inline void DofMap::enforce_constraints_on_jacobian
2375 : (const NonlinearImplicitSystem &,
2376 : SparseMatrix<Number> *) const {}
2377 :
2378 : #endif // LIBMESH_ENABLE_CONSTRAINTS
2379 :
2380 :
2381 :
2382 : inline
2383 : void DofMap::set_constrained_sparsity_construction(bool use_constraints)
2384 : {
2385 : // This got only partly finished...
2386 : if (use_constraints)
2387 : libmesh_not_implemented();
2388 :
2389 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
2390 : _constrained_sparsity_construction = use_constraints;
2391 : #endif
2392 : libmesh_ignore(use_constraints);
2393 : }
2394 :
2395 : inline
2396 : void DofMap::full_sparsity_pattern_needed()
2397 : {
2398 : need_full_sparsity_pattern = true;
2399 : }
2400 :
2401 : inline
2402 : bool DofMap::constrained_sparsity_construction()
2403 : {
2404 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
2405 : return _constrained_sparsity_construction;
2406 : #else
2407 : return true;
2408 : #endif
2409 : }
2410 :
2411 : inline
2412 : void DofMap::should_p_refine(const unsigned int g, const bool p_refine)
2413 : {
2414 : #ifdef LIBMESH_ENABLE_AMR
2415 : if (p_refine)
2416 : {
2417 : auto it = _dont_p_refine.find(g);
2418 : if (it != _dont_p_refine.end())
2419 : _dont_p_refine.erase(it);
2420 : }
2421 : else
2422 : _dont_p_refine.insert(g);
2423 : #else
2424 : libmesh_ignore(g, p_refine);
2425 : #endif
2426 : }
2427 :
2428 : inline
2429 10103642 : bool DofMap::should_p_refine(const unsigned int g) const
2430 : {
2431 : #ifdef LIBMESH_ENABLE_AMR
2432 10895767 : return !_dont_p_refine.count(g);
2433 : #else
2434 : libmesh_ignore(g);
2435 : return false;
2436 : #endif
2437 : }
2438 :
2439 : inline
2440 1181439 : unsigned int DofMap::var_group_from_var_number(const unsigned int var_num) const
2441 : {
2442 1181439 : libmesh_assert(var_num < n_variables());
2443 19152783 : return libmesh_map_find(_var_to_vg, var_num);
2444 : }
2445 :
2446 : inline
2447 6202821 : bool DofMap::should_p_refine_var(const unsigned int var) const
2448 : {
2449 : #ifdef LIBMESH_ENABLE_AMR
2450 6202821 : const auto vg = this->var_group_from_var_number(var);
2451 6202821 : return !_dont_p_refine.count(vg);
2452 : #else
2453 : libmesh_ignore(var);
2454 : return false;
2455 : #endif
2456 : }
2457 :
2458 : template <typename FieldDofsFunctor>
2459 264031877 : void DofMap::_dof_indices (const Elem & elem,
2460 : int p_level,
2461 : std::vector<dof_id_type> & di,
2462 : const unsigned int vg,
2463 : const unsigned int vig,
2464 : const Node * const * nodes,
2465 : unsigned int n_nodes,
2466 : const unsigned int v,
2467 : #ifdef DEBUG
2468 : std::size_t & tot_size,
2469 : #endif
2470 : FieldDofsFunctor field_dofs_functor) const
2471 : {
2472 264031877 : const VariableGroup & var = this->variable_group(vg);
2473 :
2474 264031877 : if (var.active_on_subdomain(elem.subdomain_id()))
2475 : {
2476 263893612 : const ElemType type = elem.type();
2477 23386464 : const unsigned int sys_num = this->sys_number();
2478 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2479 30594903 : const bool is_inf = elem.infinite();
2480 : #endif
2481 :
2482 : const bool extra_hanging_dofs =
2483 263893612 : FEInterface::extra_hanging_dofs(var.type());
2484 :
2485 263893612 : FEType fe_type = var.type();
2486 :
2487 23386464 : const bool add_p_level =
2488 : #ifdef LIBMESH_ENABLE_AMR
2489 23385196 : !_dont_p_refine.count(vg);
2490 : #else
2491 : false;
2492 : #endif
2493 :
2494 : #ifdef DEBUG
2495 : // The number of dofs per element is non-static for subdivision FE
2496 23385196 : if (var.type().family == SUBDIVISION)
2497 3936 : tot_size += n_nodes;
2498 : else
2499 : // FIXME: Is the passed-in p_level just elem.p_level()? If so,
2500 : // this seems redundant.
2501 23381260 : tot_size += FEInterface::n_dofs(fe_type, add_p_level*p_level, &elem);
2502 : #endif
2503 :
2504 : // The total Order is not required when getting the function
2505 : // pointer, it is only needed when the function is called (see
2506 : // below).
2507 : const FEInterface::n_dofs_at_node_ptr ndan =
2508 263893612 : FEInterface::n_dofs_at_node_function(fe_type, &elem);
2509 :
2510 : // Get the node-based DOF numbers
2511 1567152390 : for (unsigned int n=0; n != n_nodes; n++)
2512 : {
2513 1303258778 : const Node & node = *nodes[n];
2514 :
2515 : // Cache the intermediate lookups that are common to every
2516 : // component
2517 : #ifdef DEBUG
2518 : const std::pair<unsigned int, unsigned int>
2519 118482114 : vg_and_offset = node.var_to_vg_and_offset(sys_num,v);
2520 118482114 : libmesh_assert_equal_to (vg, vg_and_offset.first);
2521 118482114 : libmesh_assert_equal_to (vig, vg_and_offset.second);
2522 : #endif
2523 1303258778 : const unsigned int n_comp = node.n_comp_group(sys_num,vg);
2524 :
2525 : // There is a potential problem with h refinement. Imagine a
2526 : // quad9 that has a linear FE on it. Then, on the hanging side,
2527 : // it can falsely identify a DOF at the mid-edge node. This is why
2528 : // we go through FEInterface instead of node.n_comp() directly.
2529 1303258778 : const unsigned int nc =
2530 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2531 57886866 : is_inf ?
2532 82442 : FEInterface::n_dofs_at_node(fe_type, add_p_level*p_level, &elem, n) :
2533 : #endif
2534 1303185060 : ndan (type, fe_type.order + add_p_level*p_level, n);
2535 :
2536 : // If this is a non-vertex on a hanging node with extra
2537 : // degrees of freedom, we use the non-vertex dofs (which
2538 : // come in reverse order starting from the end, to
2539 : // simplify p refinement)
2540 1303258778 : if (extra_hanging_dofs && !elem.is_vertex(n))
2541 : {
2542 24870166 : const int dof_offset = n_comp - nc;
2543 :
2544 : // We should never have fewer dofs than necessary on a
2545 : // node unless we're getting indices on a parent element,
2546 : // and we should never need the indices on such a node
2547 24870166 : if (dof_offset < 0)
2548 : {
2549 0 : libmesh_assert(!elem.active());
2550 0 : di.resize(di.size() + nc, DofObject::invalid_id);
2551 : }
2552 : else
2553 56900377 : for (int i=int(n_comp)-1; i>=dof_offset; i--)
2554 : {
2555 2677598 : const dof_id_type d =
2556 32030211 : node.dof_number(sys_num, vg, vig, i, n_comp);
2557 2677598 : libmesh_assert_not_equal_to (d, DofObject::invalid_id);
2558 32030211 : field_dofs_functor(elem, n, v, di, d);
2559 : }
2560 : }
2561 : // If this is a vertex or an element without extra hanging
2562 : // dofs, our dofs come in forward order coming from the
2563 : // beginning
2564 : else
2565 : {
2566 : // We have a good component index only if it's being
2567 : // used on this FE type (nc) *and* it's available on
2568 : // this DofObject (n_comp).
2569 1278388612 : const unsigned int good_nc = std::min(n_comp, nc);
2570 2474735368 : for (unsigned int i=0; i!=good_nc; ++i)
2571 : {
2572 106364690 : const dof_id_type d =
2573 1196341872 : node.dof_number(sys_num, vg, vig, i, n_comp);
2574 106364690 : libmesh_assert_not_equal_to (d, DofObject::invalid_id);
2575 106364690 : libmesh_assert_less (d, this->n_dofs());
2576 1196346756 : field_dofs_functor(elem, n, v, di, d);
2577 : }
2578 :
2579 : // With fewer good component indices than we need, e.g.
2580 : // due to subdomain expansion, the remaining expected
2581 : // indices are marked invalid.
2582 1278388612 : if (n_comp < nc)
2583 0 : for (unsigned int i=n_comp; i!=nc; ++i)
2584 0 : di.push_back(DofObject::invalid_id);
2585 : }
2586 : }
2587 :
2588 : // If there are any element-based DOF numbers, get them
2589 263893612 : const unsigned int nc = FEInterface::n_dofs_per_elem(fe_type, add_p_level*p_level, &elem);
2590 :
2591 : // We should never have fewer dofs than necessary on an
2592 : // element unless we're getting indices on a parent element
2593 : // (and we should never need those indices) or off-domain for a
2594 : // subdomain-restricted variable (where invalid_id is the
2595 : // correct thing to return)
2596 263893612 : if (nc != 0)
2597 : {
2598 14798645 : const unsigned int n_comp = elem.n_comp_group(sys_num,vg);
2599 14798645 : if (elem.n_systems() > sys_num && nc <= n_comp)
2600 : {
2601 90064574 : for (unsigned int i=0; i<nc; i++)
2602 : {
2603 6779581 : const dof_id_type d =
2604 75263241 : elem.dof_number(sys_num, vg, vig, i, n_comp);
2605 6779581 : libmesh_assert_not_equal_to (d, DofObject::invalid_id);
2606 :
2607 75265929 : field_dofs_functor(elem, invalid_uint, v, di, d);
2608 : }
2609 : }
2610 : else
2611 : {
2612 0 : libmesh_assert(!elem.active() || fe_type.family == LAGRANGE || fe_type.family == SUBDIVISION);
2613 0 : di.resize(di.size() + nc, DofObject::invalid_id);
2614 : }
2615 : }
2616 : }
2617 264031877 : }
2618 :
2619 :
2620 :
2621 : template <typename ScalarDofsFunctor, typename FieldDofsFunctor>
2622 110983826 : void DofMap::dof_indices (const Elem * const elem,
2623 : std::vector<dof_id_type> & di,
2624 : const unsigned int vn,
2625 : ScalarDofsFunctor scalar_dofs_functor,
2626 : FieldDofsFunctor field_dofs_functor,
2627 : int p_level) const
2628 : {
2629 : // We now allow elem==nullptr to request just SCALAR dofs
2630 : // libmesh_assert(elem);
2631 :
2632 10723403 : LOG_SCOPE("dof_indices()", "DofMap");
2633 :
2634 : // Clear the DOF indices vector
2635 10723403 : di.clear();
2636 :
2637 : // Use the default p refinement level?
2638 110983826 : if (p_level == -12345)
2639 110358818 : p_level = elem ? elem->p_level() : 0;
2640 :
2641 110985126 : const unsigned int vg = this->_variable_group_numbers[vn];
2642 10723403 : const VariableGroup & var = this->variable_group(vg);
2643 110983826 : const unsigned int vig = vn - var.number();
2644 :
2645 : #ifdef DEBUG
2646 : // Check that sizes match in DEBUG mode
2647 10723403 : std::size_t tot_size = 0;
2648 : #endif
2649 :
2650 110983826 : if (elem && elem->type() == TRI3SUBDIVISION)
2651 : {
2652 : // Subdivision surface FE require the 1-ring around elem
2653 2712 : const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
2654 :
2655 : // Ghost subdivision elements have no real dofs
2656 29832 : if (!sd_elem->is_ghost())
2657 : {
2658 : // Determine the nodes contributing to element elem
2659 4608 : std::vector<const Node *> elem_nodes;
2660 25344 : MeshTools::Subdivision::find_one_ring(sd_elem, elem_nodes);
2661 :
2662 25344 : _dof_indices(*elem, p_level, di, vg, vig, elem_nodes.data(),
2663 : cast_int<unsigned int>(elem_nodes.size()), vn,
2664 : #ifdef DEBUG
2665 : tot_size,
2666 : #endif
2667 : field_dofs_functor);
2668 : }
2669 :
2670 29832 : return;
2671 : }
2672 :
2673 : // Get the dof numbers
2674 111051945 : if (var.type().family == SCALAR &&
2675 1027177 : (!elem ||
2676 1027209 : var.active_on_subdomain(elem->subdomain_id())))
2677 : {
2678 : #ifdef DEBUG
2679 97951 : tot_size += var.type().order;
2680 : #endif
2681 195902 : std::vector<dof_id_type> di_new;
2682 1027209 : this->SCALAR_dof_indices(di_new,vn);
2683 1027209 : scalar_dofs_functor(*elem, di, di_new);
2684 : }
2685 109926785 : else if (elem)
2686 109926785 : _dof_indices(*elem, p_level, di, vg, vig, elem->get_nodes(),
2687 109926785 : elem->n_nodes(), vn,
2688 : #ifdef DEBUG
2689 : tot_size,
2690 : #endif
2691 : field_dofs_functor);
2692 :
2693 : #ifdef DEBUG
2694 10720691 : libmesh_assert_equal_to (tot_size, di.size());
2695 : #endif
2696 : }
2697 :
2698 : inline
2699 10 : StaticCondensationDofMap & DofMap::get_static_condensation()
2700 : {
2701 10 : libmesh_assert(_sc);
2702 10 : return *_sc;
2703 : }
2704 :
2705 : } // namespace libMesh
2706 :
2707 : #endif // LIBMESH_DOF_MAP_H
|
