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_SYSTEM_H
21 : #define LIBMESH_SYSTEM_H
22 :
23 : // Local Includes
24 : #include "libmesh/elem_range.h"
25 : #include "libmesh/enum_subset_solve_mode.h" // SUBSET_ZERO
26 : #include "libmesh/enum_parallel_type.h" // PARALLEL
27 : #include "libmesh/fem_function_base.h"
28 : #include "libmesh/libmesh_common.h"
29 : #include "libmesh/parallel_object.h"
30 : #include "libmesh/parameters.h"
31 : #include "libmesh/qoi_set.h"
32 : #include "libmesh/reference_counted_object.h"
33 : #include "libmesh/tensor_value.h" // For point_hessian
34 : #include "libmesh/enum_matrix_build_type.h" // AUTOMATIC
35 :
36 : // C++ includes
37 : #include <cstddef>
38 : #include <set>
39 : #include <string>
40 : #include <string_view>
41 : #include <vector>
42 : #include <memory>
43 :
44 : // This define may be useful in --disable-optional builds when it is
45 : // possible that libmesh will not have any solvers available.
46 : #if defined(LIBMESH_HAVE_PETSC) || \
47 : defined(LIBMESH_HAVE_EIGEN) || \
48 : defined(LIBMESH_HAVE_LASPACK) || \
49 : defined(LIBMESH_TRILINOS_HAVE_AZTECOO)
50 : #define LIBMESH_HAVE_SOLVER 1
51 : #endif
52 :
53 : namespace libMesh
54 : {
55 :
56 : // Forward Declarations
57 : class System;
58 : class EquationSystems;
59 : class MeshBase;
60 : class Xdr;
61 : class DofMap;
62 : template <typename Output> class FunctionBase;
63 : class ParameterVector;
64 : class Point;
65 : class SensitivityData;
66 : class Matrix;
67 : template <typename T> class NumericVector;
68 : template <typename T> class SparseMatrix;
69 : template <typename T> class VectorValue;
70 : typedef VectorValue<Number> NumberVectorValue;
71 : typedef NumberVectorValue Gradient;
72 : class SystemSubset;
73 : class FEType;
74 : class SystemNorm;
75 : enum FEMNormType : int;
76 : class Variable;
77 : class VariableGroup;
78 :
79 : /**
80 : * \brief Manages consistently variables, degrees of freedom, and coefficient
81 : * vectors.
82 : *
83 : * This is the base class for classes which contain information related to any
84 : * physical process that might be simulated. Such information may range from
85 : * the actual solution values to algorithmic flags that may be used to control
86 : * the numerical methods employed. In general, use an EquationSystems
87 : * object to handle one or more of the children of this class.
88 : *
89 : * Especially, this class manages the variables of the differential equation,
90 : * the coefficient vectors and the \p DofMap, and ensures that these are
91 : * consistent. It provides storage for the solution. Furthermore, (de-)
92 : * serialization functionality is provided.
93 : *
94 : * \author Benjamin S. Kirk
95 : * \date 2003-2004
96 : */
97 : class System : public ReferenceCountedObject<System>,
98 : public ParallelObject
99 : {
100 : public:
101 :
102 : /**
103 : * Constructor. Optionally initializes required
104 : * data structures.
105 : */
106 : System (EquationSystems & es,
107 : const std::string & name,
108 : const unsigned int number);
109 :
110 : /**
111 : * Special functions.
112 : * - The copy constructor and assignment operator were previously
113 : * declared as private libmesh_not_implemented() functions, this
114 : * is the C++11 way of achieving the same effect.
115 : * - The System holds references to Mesh and EquationSystems
116 : * objects, therefore it can't be default move-assigned.
117 : * - This class _can_ be default move constructed.
118 : * - The destructor dies with an error in dbg/devel modes if libMesh::closed()
119 : */
120 : System (const System &) = delete;
121 : System & operator= (const System &) = delete;
122 : System (System &&) = default;
123 : System & operator= (System &&) = delete;
124 : virtual ~System ();
125 :
126 : /**
127 : * Abstract base class to be used for system initialization.
128 : * A user class derived from this class may be used to
129 : * initialize the system values by attaching an object
130 : * with the method \p attach_init_object.
131 : */
132 : class Initialization
133 : {
134 : public:
135 : /**
136 : * Destructor. Virtual because we will have virtual functions.
137 : */
138 : virtual ~Initialization () = default;
139 :
140 : /**
141 : * Initialization function. This function will be called
142 : * to initialize the system values upon creation and must
143 : * be provided by the user in a derived class.
144 : */
145 : virtual void initialize () = 0;
146 : };
147 :
148 :
149 :
150 : /**
151 : * Abstract base class to be used for system assembly.
152 : * A user class derived from this class may be used to
153 : * assemble the system by attaching an object
154 : * with the method \p attach_assemble_object.
155 : */
156 : class Assembly
157 : {
158 : public:
159 : /**
160 : * Destructor. Virtual because we will have virtual functions.
161 : */
162 : virtual ~Assembly () = default;
163 :
164 : /**
165 : * Assembly function. This function will be called
166 : * to assemble the system prior to a solve and must
167 : * be provided by the user in a derived class.
168 : */
169 : virtual void assemble () = 0;
170 : };
171 :
172 :
173 :
174 : /**
175 : * Abstract base class to be used for system constraints.
176 : * A user class derived from this class may be used to
177 : * constrain the system by attaching an object
178 : * with the method \p attach_constraint_object.
179 : */
180 : class Constraint
181 : {
182 : public:
183 : /**
184 : * Destructor. Virtual because we will have virtual functions.
185 : */
186 0 : virtual ~Constraint () = default;
187 :
188 : /**
189 : * Constraint function. This function will be called
190 : * to constrain the system prior to a solve and must
191 : * be provided by the user in a derived class.
192 : */
193 : virtual void constrain () = 0;
194 : };
195 :
196 :
197 :
198 : /**
199 : * Abstract base class to be used for quantities of interest.
200 : * A user class derived from this class may be used to
201 : * compute quantities of interest by attaching an object
202 : * with the method \p attach_QOI_object.
203 : */
204 : class QOI
205 : {
206 : public:
207 : /**
208 : * Destructor. Virtual because we will have virtual functions.
209 : */
210 : virtual ~QOI () = default;
211 :
212 : /**
213 : * Quantity of interest function. This function will be called
214 : * to compute quantities of interest and must be provided by the
215 : * user in a derived class.
216 : */
217 : virtual void qoi (const QoISet & qoi_indices) = 0;
218 : };
219 :
220 :
221 :
222 : /**
223 : * Abstract base class to be used for derivatives of quantities
224 : * of interest. A user class derived from this class may be used
225 : * to compute quantities of interest by attaching an object
226 : * with the method \p attach_QOI_derivative_object.
227 : */
228 : class QOIDerivative
229 : {
230 : public:
231 : /**
232 : * Destructor. Virtual because we will have virtual functions.
233 : */
234 : virtual ~QOIDerivative () = default;
235 :
236 : /**
237 : * Quantity of interest derivative function. This function will
238 : * be called to compute derivatives of quantities of interest and
239 : * must be provided by the user in a derived class.
240 : */
241 : virtual void qoi_derivative (const QoISet & qoi_indices,
242 : bool include_liftfunc,
243 : bool apply_constraints) = 0;
244 : };
245 :
246 : /**
247 : * The type of system.
248 : */
249 : typedef System sys_type;
250 :
251 : /**
252 : * \returns A reference to *this.
253 : */
254 : sys_type & system () { return *this; }
255 :
256 : /**
257 : * Clear all the data structures associated with
258 : * the system.
259 : */
260 : virtual void clear ();
261 :
262 : /**
263 : * Initializes degrees of freedom on the current mesh.
264 : * Sets the
265 : */
266 : void init ();
267 :
268 : /**
269 : * Reinitializes degrees of freedom and other required data on the
270 : * current mesh.
271 : *
272 : * \note The matrix is not initialized at this time since it may not
273 : * be required for all applications. Should be overridden in derived
274 : * classes.
275 : */
276 : virtual void reinit ();
277 :
278 : /**
279 : * Reinitializes the constraints for this system.
280 : *
281 : * Also prepares the send_list, whether or not constraints have
282 : * changed.
283 : */
284 : virtual void reinit_constraints ();
285 :
286 : /**
287 : * Reinitializes the system with a new mesh.
288 : */
289 : virtual void reinit_mesh();
290 :
291 : /**
292 : * \returns \p true iff this system has been initialized.
293 : */
294 : bool is_initialized() const;
295 :
296 : /**
297 : * Update the local values to reflect the solution
298 : * on neighboring processors.
299 : */
300 : virtual void update ();
301 :
302 : /**
303 : * Prepares \p matrix and \p _dof_map for matrix assembly.
304 : * Does not actually assemble anything. For matrix assembly,
305 : * use the \p assemble() in derived classes.
306 : * Should be overridden in derived classes.
307 : */
308 : virtual void assemble ();
309 :
310 : /**
311 : * Calls user qoi function.
312 : * Can be overridden in derived classes.
313 : */
314 : virtual void assemble_qoi (const QoISet & qoi_indices = QoISet());
315 :
316 : /**
317 : * Calls user qoi derivative function.
318 : * Can be overridden in derived classes.
319 : */
320 : virtual void assemble_qoi_derivative (const QoISet & qoi_indices = QoISet(),
321 : bool include_liftfunc = true,
322 : bool apply_constraints = true);
323 :
324 : /**
325 : * Calls residual parameter derivative function.
326 : *
327 : * Library subclasses use finite differences by default.
328 : *
329 : * This should assemble the sensitivity rhs vectors to hold
330 : * -(partial R / partial p_i), making them ready to solve
331 : * the forward sensitivity equation.
332 : *
333 : * This method is only implemented in some derived classes.
334 : */
335 : virtual void assemble_residual_derivatives (const ParameterVector & parameters);
336 :
337 : /**
338 : * After calling this method, any solve will be restricted to the
339 : * given subdomain. To disable this mode, call this method with \p
340 : * subset being a \p nullptr.
341 : */
342 : virtual void restrict_solve_to (const SystemSubset * subset,
343 : const SubsetSolveMode subset_solve_mode=SUBSET_ZERO);
344 :
345 : /**
346 : * Solves the system. Should be overridden in derived systems.
347 : */
348 0 : virtual void solve () {}
349 :
350 : /**
351 : * Solves the sensitivity system, for the provided parameters.
352 : * Must be overridden in derived systems.
353 : *
354 : * \returns A pair with the total number of linear iterations
355 : * performed and the (sum of the) final residual norms
356 : *
357 : * This method is only implemented in some derived classes.
358 : */
359 : virtual std::pair<unsigned int, Real>
360 : sensitivity_solve (const ParameterVector & parameters);
361 :
362 : /**
363 : * Assembles & solves the linear system(s) (dR/du)*u_w = sum(w_p*-dR/dp), for
364 : * those parameters p contained within \p parameters weighted by the
365 : * values w_p found within \p weights.
366 : *
367 : * \returns A pair with the total number of linear iterations
368 : * performed and the (sum of the) final residual norms
369 : *
370 : * This method is only implemented in some derived classes.
371 : */
372 : virtual std::pair<unsigned int, Real>
373 : weighted_sensitivity_solve (const ParameterVector & parameters,
374 : const ParameterVector & weights);
375 :
376 : /**
377 : * Solves the adjoint system, for the specified qoi indices, or for
378 : * every qoi if \p qoi_indices is nullptr. Must be overridden in
379 : * derived systems.
380 : *
381 : * \returns A pair with the total number of linear iterations
382 : * performed and the (sum of the) final residual norms
383 : *
384 : * This method is only implemented in some derived classes.
385 : */
386 : virtual std::pair<unsigned int, Real>
387 : adjoint_solve (const QoISet & qoi_indices = QoISet());
388 :
389 : /**
390 : * Assembles & solves the linear system(s)
391 : * (dR/du)^T*z_w = sum(w_p*(d^2q/dudp - d^2R/dudp*z)), for those
392 : * parameters p contained within \p parameters, weighted by the
393 : * values w_p found within \p weights.
394 : *
395 : * Assumes that adjoint_solve has already calculated z for each qoi
396 : * in \p qoi_indices.
397 : *
398 : * \returns A pair with the total number of linear iterations
399 : * performed and the (sum of the) final residual norms
400 : *
401 : * This method is only implemented in some derived classes.
402 : */
403 : virtual std::pair<unsigned int, Real>
404 : weighted_sensitivity_adjoint_solve (const ParameterVector & parameters,
405 : const ParameterVector & weights,
406 : const QoISet & qoi_indices = QoISet());
407 : /**
408 : * Accessor for the adjoint_already_solved boolean
409 : */
410 784 : bool is_adjoint_already_solved() const
411 25384 : { return adjoint_already_solved;}
412 :
413 : /**
414 : * Setter for the adjoint_already_solved boolean
415 : */
416 : void set_adjoint_already_solved(bool setting)
417 : { adjoint_already_solved = setting;}
418 :
419 :
420 : /**
421 : * Solves for the derivative of each of the system's quantities of
422 : * interest q in \p qoi[qoi_indices] with respect to each parameter
423 : * in \p parameters, placing the result for qoi \p i and parameter
424 : * \p j into \p sensitivities[i][j].
425 : *
426 : * \note \p parameters is a const vector, not a vector-of-const;
427 : * parameter values in this vector need to be mutable for finite
428 : * differencing to work.
429 : *
430 : * Automatically chooses the forward method for problems with more
431 : * quantities of interest than parameters, or the adjoint method
432 : * otherwise.
433 : *
434 : * This method is only usable in derived classes which override
435 : * an implementation.
436 : */
437 : virtual void qoi_parameter_sensitivity (const QoISet & qoi_indices,
438 : const ParameterVector & parameters,
439 : SensitivityData & sensitivities);
440 :
441 : /**
442 : * Solves for parameter sensitivities using the adjoint method.
443 : *
444 : * This method is only implemented in some derived classes.
445 : */
446 : virtual void adjoint_qoi_parameter_sensitivity (const QoISet & qoi_indices,
447 : const ParameterVector & parameters,
448 : SensitivityData & sensitivities);
449 :
450 : /**
451 : * Solves for parameter sensitivities using the forward method.
452 : *
453 : * This method is only implemented in some derived classes.
454 : */
455 : virtual void forward_qoi_parameter_sensitivity (const QoISet & qoi_indices,
456 : const ParameterVector & parameters,
457 : SensitivityData & sensitivities);
458 :
459 : /**
460 : * For each of the system's quantities of interest q in
461 : * \p qoi[qoi_indices], and for a vector of parameters p, the
462 : * parameter sensitivity Hessian H_ij is defined as
463 : * H_ij = (d^2 q)/(d p_i d p_j)
464 : * This Hessian is the output of this method, where for each q_i,
465 : * H_jk is stored in \p hessian.second_derivative(i,j,k).
466 : *
467 : * This method is only implemented in some derived classes.
468 : */
469 : virtual void qoi_parameter_hessian(const QoISet & qoi_indices,
470 : const ParameterVector & parameters,
471 : SensitivityData & hessian);
472 :
473 : /**
474 : * For each of the system's quantities of interest q in
475 : * \p qoi[qoi_indices], and for a vector of parameters p, the
476 : * parameter sensitivity Hessian H_ij is defined as
477 : * H_ij = (d^2 q)/(d p_i d p_j)
478 : * The Hessian-vector product, for a vector v_k in parameter space, is
479 : * S_j = H_jk v_k
480 : * This product is the output of this method, where for each q_i,
481 : * S_j is stored in \p sensitivities[i][j].
482 : *
483 : * This method is only implemented in some derived classes.
484 : */
485 : virtual void qoi_parameter_hessian_vector_product(const QoISet & qoi_indices,
486 : const ParameterVector & parameters,
487 : const ParameterVector & vector,
488 : SensitivityData & product);
489 :
490 : /**
491 : * \returns \p true when the other system contains
492 : * identical data, up to the given threshold. Outputs
493 : * some diagnostic info when \p verbose is set.
494 : */
495 : virtual bool compare (const System & other_system,
496 : const Real threshold,
497 : const bool verbose) const;
498 :
499 : /**
500 : * \returns The system name.
501 : */
502 : const std::string & name () const;
503 :
504 : /**
505 : * \returns The type of system, helpful in identifying
506 : * which system type to use when reading equation system
507 : * data from file. Should be overridden in derived classes.
508 : */
509 502 : virtual std::string system_type () const { return "Basic"; }
510 :
511 : /**
512 : * Projects arbitrary functions onto the current solution.
513 : * The function value \p f and its gradient \p g are
514 : * user-provided cloneable functors.
515 : * A gradient \p g is only required/used for projecting onto finite
516 : * element spaces with continuous derivatives.
517 : */
518 : void project_solution (FunctionBase<Number> * f,
519 : FunctionBase<Gradient> * g = nullptr) const;
520 :
521 : /**
522 : * Projects arbitrary functions onto the current solution.
523 : * The function value \p f and its gradient \p g are
524 : * user-provided cloneable functors.
525 : * A gradient \p g is only required/used for projecting onto finite
526 : * element spaces with continuous derivatives.
527 : */
528 : void project_solution (FEMFunctionBase<Number> * f,
529 : FEMFunctionBase<Gradient> * g = nullptr) const;
530 :
531 : /**
532 : * Projects arbitrary functions onto the current solution.
533 : * The function value \p fptr and its gradient \p gptr are
534 : * represented by function pointers.
535 : * A gradient \p gptr is only required/used for projecting onto
536 : * finite element spaces with continuous derivatives.
537 : */
538 : typedef Number (*ValueFunctionPointer)(const Point & p,
539 : const Parameters & Parameters,
540 : const std::string & sys_name,
541 : const std::string & unknown_name);
542 : typedef Gradient (*GradientFunctionPointer)(const Point & p,
543 : const Parameters & parameters,
544 : const std::string & sys_name,
545 : const std::string & unknown_name);
546 : void project_solution (ValueFunctionPointer fptr,
547 : GradientFunctionPointer gptr,
548 : const Parameters & parameters) const;
549 :
550 : /**
551 : * Projects arbitrary functions onto a vector of degree of freedom
552 : * values for the current system.
553 : * The function value \p f and its gradient \p g are
554 : * user-provided cloneable functors.
555 : * A gradient \p g is only required/used for projecting onto finite
556 : * element spaces with continuous derivatives.
557 : *
558 : * Constrain the new vector using the requested adjoint rather than
559 : * primal constraints if is_adjoint is non-negative.
560 : */
561 : void project_vector (NumericVector<Number> & new_vector,
562 : FunctionBase<Number> * f,
563 : FunctionBase<Gradient> * g = nullptr,
564 : int is_adjoint = -1) const;
565 :
566 : /**
567 : * Projects arbitrary functions onto a vector of degree of freedom
568 : * values for the current system.
569 : * The function value \p f and its gradient \p g are
570 : * user-provided cloneable functors.
571 : * A gradient \p g is only required/used for projecting onto finite
572 : * element spaces with continuous derivatives.
573 : *
574 : * Constrain the new vector using the requested adjoint rather than
575 : * primal constraints if is_adjoint is non-negative.
576 : */
577 : void project_vector (NumericVector<Number> & new_vector,
578 : FEMFunctionBase<Number> * f,
579 : FEMFunctionBase<Gradient> * g = nullptr,
580 : int is_adjoint = -1) const;
581 :
582 : /**
583 : * Projects arbitrary functions onto a vector of degree of freedom
584 : * values for the current system.
585 : * The function value \p fptr and its gradient \p gptr are
586 : * represented by function pointers.
587 : * A gradient \p gptr is only required/used for projecting onto
588 : * finite element spaces with continuous derivatives.
589 : *
590 : * Constrain the new vector using the requested adjoint rather than
591 : * primal constraints if is_adjoint is non-negative.
592 : */
593 : void project_vector (ValueFunctionPointer fptr,
594 : GradientFunctionPointer gptr,
595 : const Parameters & parameters,
596 : NumericVector<Number> & new_vector,
597 : int is_adjoint = -1) const;
598 :
599 : /**
600 : * Projects arbitrary boundary functions onto a vector of degree of
601 : * freedom values for the current system.
602 : * Only degrees of freedom which affect the function's trace on a
603 : * boundary in the set \p b are affected.
604 : * Only degrees of freedom associated with the variables listed in
605 : * the vector \p variables are projected.
606 : * The function value \p f and its gradient \p g are
607 : * user-provided cloneable functors.
608 : * A gradient \p g is only required/used for projecting onto finite
609 : * element spaces with continuous derivatives.
610 : */
611 : void boundary_project_solution (const std::set<boundary_id_type> & b,
612 : const std::vector<unsigned int> & variables,
613 : FunctionBase<Number> * f,
614 : FunctionBase<Gradient> * g = nullptr);
615 :
616 : /**
617 : * Projects arbitrary boundary functions onto a vector of degree of
618 : * freedom values for the current system.
619 : * Only degrees of freedom which affect the function's trace on a
620 : * boundary in the set \p b are affected.
621 : * Only degrees of freedom associated with the variables listed in
622 : * the vector \p variables are projected.
623 : * The function value \p fptr and its gradient \p gptr are
624 : * represented by function pointers.
625 : * A gradient \p gptr is only required/used for projecting onto
626 : * finite element spaces with continuous derivatives.
627 : */
628 : void boundary_project_solution (const std::set<boundary_id_type> & b,
629 : const std::vector<unsigned int> & variables,
630 : ValueFunctionPointer fptr,
631 : GradientFunctionPointer gptr,
632 : const Parameters & parameters);
633 :
634 : /**
635 : * Projects arbitrary boundary functions onto a vector of degree of
636 : * freedom values for the current system.
637 : * Only degrees of freedom which affect the function's trace on a
638 : * boundary in the set \p b are affected.
639 : * Only degrees of freedom associated with the variables listed in
640 : * the vector \p variables are projected.
641 : * The function value \p f and its gradient \p g are
642 : * user-provided cloneable functors.
643 : * A gradient \p g is only required/used for projecting onto finite
644 : * element spaces with continuous derivatives.
645 : *
646 : * Constrain the new vector using the requested adjoint rather than
647 : * primal constraints if is_adjoint is non-negative.
648 : */
649 : void boundary_project_vector (const std::set<boundary_id_type> & b,
650 : const std::vector<unsigned int> & variables,
651 : NumericVector<Number> & new_vector,
652 : FunctionBase<Number> * f,
653 : FunctionBase<Gradient> * g = nullptr,
654 : int is_adjoint = -1) const;
655 :
656 : /**
657 : * Projects arbitrary boundary functions onto a vector of degree of
658 : * freedom values for the current system.
659 : * Only degrees of freedom which affect the function's trace on a
660 : * boundary in the set \p b are affected.
661 : * Only degrees of freedom associated with the variables listed in
662 : * the vector \p variables are projected.
663 : * The function value \p fptr and its gradient \p gptr are
664 : * represented by function pointers.
665 : * A gradient \p gptr is only required/used for projecting onto
666 : * finite element spaces with continuous derivatives.
667 : *
668 : * Constrain the new vector using the requested adjoint rather than
669 : * primal constraints if is_adjoint is non-negative.
670 : */
671 : void boundary_project_vector (const std::set<boundary_id_type> & b,
672 : const std::vector<unsigned int> & variables,
673 : ValueFunctionPointer fptr,
674 : GradientFunctionPointer gptr,
675 : const Parameters & parameters,
676 : NumericVector<Number> & new_vector,
677 : int is_adjoint = -1) const;
678 :
679 : /**
680 : * \returns The system number.
681 : */
682 : unsigned int number () const;
683 :
684 : /**
685 : * Fill the input vector \p global_soln so that it contains
686 : * the global solution on all processors.
687 : * Requires communication with all other processors.
688 : */
689 : void update_global_solution (std::vector<Number> & global_soln) const;
690 :
691 : /**
692 : * Fill the input vector \p global_soln so that it contains
693 : * the global solution on processor \p dest_proc.
694 : * Requires communication with all other processors.
695 : */
696 : void update_global_solution (std::vector<Number> & global_soln,
697 : const processor_id_type dest_proc) const;
698 :
699 : /**
700 : * \returns A constant reference to this systems's \p _mesh.
701 : */
702 : const MeshBase & get_mesh() const;
703 :
704 : /**
705 : * \returns A reference to this systems's \p _mesh.
706 : */
707 : MeshBase & get_mesh();
708 :
709 : /**
710 : * \returns A constant reference to this system's \p _dof_map.
711 : */
712 : const DofMap & get_dof_map() const;
713 :
714 : /**
715 : * \returns A writable reference to this system's \p _dof_map.
716 : */
717 : DofMap & get_dof_map();
718 :
719 : /**
720 : * \returns A constant reference to this system's parent EquationSystems object.
721 : */
722 558896 : const EquationSystems & get_equation_systems() const { return _equation_systems; }
723 :
724 : /**
725 : * \returns A reference to this system's parent EquationSystems object.
726 : */
727 219069 : EquationSystems & get_equation_systems() { return _equation_systems; }
728 :
729 : /**
730 : * \returns \p true if the system is active, \p false otherwise.
731 : * An active system will be solved.
732 : */
733 : bool active () const;
734 :
735 : /**
736 : * Activates the system. Only active systems are solved.
737 : */
738 : void activate ();
739 :
740 : /**
741 : * Deactivates the system. Only active systems are solved.
742 : */
743 : void deactivate ();
744 :
745 : /**
746 : * Sets the system to be "basic only": i.e. advanced system
747 : * components such as ImplicitSystem matrices may not be
748 : * initialized. This is useful for efficiency in certain utility
749 : * programs that never use System::solve(). This method must be
750 : * called after the System or derived class is created but before it
751 : * is initialized; e.g. from within EquationSystems::read()
752 : */
753 : void set_basic_system_only ();
754 :
755 : /**
756 : * Vector iterator typedefs.
757 : */
758 : typedef std::map<std::string, std::unique_ptr<NumericVector<Number>>, std::less<>>::iterator vectors_iterator;
759 : typedef std::map<std::string, std::unique_ptr<NumericVector<Number>>, std::less<>>::const_iterator const_vectors_iterator;
760 :
761 : /**
762 : * Beginning of vectors container
763 : */
764 : vectors_iterator vectors_begin ();
765 :
766 : /**
767 : * Beginning of vectors container
768 : */
769 : const_vectors_iterator vectors_begin () const;
770 :
771 : /**
772 : * End of vectors container
773 : */
774 : vectors_iterator vectors_end ();
775 :
776 : /**
777 : * End of vectors container
778 : */
779 : const_vectors_iterator vectors_end () const;
780 :
781 : /**
782 : * Matrix iterator typedefs.
783 : */
784 : typedef std::map<std::string, std::unique_ptr<SparseMatrix<Number>>, std::less<>>::iterator matrices_iterator;
785 : typedef std::map<std::string, std::unique_ptr<SparseMatrix<Number>>, std::less<>>::const_iterator const_matrices_iterator;
786 :
787 : /**
788 : * Beginning of matrices container
789 : */
790 : matrices_iterator matrices_begin ();
791 :
792 : /**
793 : * Beginning of matrices container
794 : */
795 : const_matrices_iterator matrices_begin () const;
796 :
797 : /**
798 : * End of matrices container
799 : */
800 : matrices_iterator matrices_end ();
801 :
802 : /**
803 : * End of matrices container
804 : */
805 : const_matrices_iterator matrices_end () const;
806 :
807 : /**
808 : * Adds the additional vector \p vec_name to this system. All the
809 : * additional vectors are similarly distributed, like the \p
810 : * solution, and initialized to zero.
811 : *
812 : * By default vectors added by add_vector are projected to changed grids by
813 : * reinit(). To zero them instead (more efficient), pass "false" as the
814 : * second argument
815 : *
816 : * If the vector already exists, the existing vector is returned.
817 : * after any upgrade to the \p projections or \p type has been made.
818 : * We *only* handle upgrades (projections false->true, or type
819 : * PARALLEL->GHOSTED) in this fashion, not downgrades, on the theory
820 : * that if two codes have differing needs we want to support the
821 : * union of those needs, not the intersection. Downgrades can only
822 : * be accomplished manually, via \p set_vector_preservation() or by
823 : * setting a vector \p type() and re-initializing.
824 : */
825 : NumericVector<Number> & add_vector (std::string_view vec_name,
826 : const bool projections=true,
827 : const ParallelType type = PARALLEL);
828 :
829 : /**
830 : * Removes the additional vector \p vec_name from this system
831 : */
832 : void remove_vector(std::string_view vec_name);
833 :
834 : /**
835 : * Tells the System whether or not to project the solution vector onto new
836 : * grids when the system is reinitialized. The solution will be projected
837 : * unless project_solution_on_reinit() = false is called.
838 : */
839 2392 : bool & project_solution_on_reinit (void)
840 2392 : { return _solution_projection; }
841 :
842 : /**
843 : * \returns \p true if this \p System has a vector associated with the
844 : * given name, \p false otherwise.
845 : */
846 : bool have_vector (std::string_view vec_name) const;
847 :
848 : /**
849 : * \returns A const pointer to the vector if this \p System has a
850 : * vector associated with the given name, \p nullptr otherwise.
851 : */
852 : const NumericVector<Number> * request_vector (std::string_view vec_name) const;
853 :
854 : /**
855 : * \returns A pointer to the vector if this \p System has a
856 : * vector associated with the given name, \p nullptr otherwise.
857 : */
858 : NumericVector<Number> * request_vector (std::string_view vec_name);
859 :
860 : /**
861 : * \returns A const pointer to this system's additional vector
862 : * number \p vec_num (where the vectors are counted starting with
863 : * 0), or \p nullptr if the system has no such vector.
864 : */
865 : const NumericVector<Number> * request_vector (const unsigned int vec_num) const;
866 :
867 : /**
868 : * \returns A writable pointer to this system's additional
869 : * vector number \p vec_num (where the vectors are counted starting
870 : * with 0), or \p nullptr if the system has no such vector.
871 : */
872 : NumericVector<Number> * request_vector (const unsigned int vec_num);
873 :
874 : /**
875 : * \returns A const reference to this system's additional vector
876 : * named \p vec_name. Access is only granted when the vector is already
877 : * properly initialized.
878 : */
879 : const NumericVector<Number> & get_vector (std::string_view vec_name) const;
880 :
881 : /**
882 : * \returns A writable reference to this system's additional vector
883 : * named \p vec_name. Access is only granted when the vector is already
884 : * properly initialized.
885 : */
886 : NumericVector<Number> & get_vector (std::string_view vec_name);
887 :
888 : /**
889 : * \returns A const reference to this system's additional vector
890 : * number \p vec_num (where the vectors are counted starting with
891 : * 0).
892 : */
893 : const NumericVector<Number> & get_vector (const unsigned int vec_num) const;
894 :
895 : /**
896 : * \returns A writable reference to this system's additional
897 : * vector number \p vec_num (where the vectors are counted starting
898 : * with 0).
899 : */
900 : NumericVector<Number> & get_vector (const unsigned int vec_num);
901 :
902 : /**
903 : * \returns The name of this system's additional vector number \p
904 : * vec_num (where the vectors are counted starting with 0).
905 : */
906 : const std::string & vector_name (const unsigned int vec_num) const;
907 :
908 : /**
909 : * \returns The name of a system vector, given a reference to that vector
910 : */
911 : const std::string & vector_name (const NumericVector<Number> & vec_reference) const;
912 :
913 : /**
914 : * Allows one to set the QoI index controlling whether the vector
915 : * identified by vec_name represents a solution from the adjoint
916 : * (qoi_num >= 0) or primal (qoi_num == -1) space. This becomes
917 : * significant if those spaces have differing heterogeneous
918 : * Dirichlet constraints.
919 : *
920 : * qoi_num == -2 can be used to indicate a vector which should not
921 : * be affected by constraints during projection operations.
922 : */
923 : void set_vector_as_adjoint (const std::string & vec_name, int qoi_num);
924 :
925 : /**
926 : * \returns The integer describing whether the vector identified by
927 : * vec_name represents a solution from an adjoint (non-negative) or
928 : * the primal (-1) space.
929 : */
930 : int vector_is_adjoint (std::string_view vec_name) const;
931 :
932 : /**
933 : * Allows one to set the boolean controlling whether the vector
934 : * identified by vec_name should be "preserved": projected to new
935 : * meshes, saved, etc.
936 : */
937 : void set_vector_preservation (const std::string & vec_name, bool preserve);
938 :
939 : /**
940 : * \returns The boolean describing whether the vector identified by
941 : * vec_name should be "preserved": projected to new meshes, saved,
942 : * etc.
943 : */
944 : bool vector_preservation (std::string_view vec_name) const;
945 :
946 : /**
947 : * \returns A reference to one of the system's adjoint solution
948 : * vectors, by default the one corresponding to the first qoi.
949 : * Creates the vector if it doesn't already exist.
950 : */
951 : NumericVector<Number> & add_adjoint_solution(unsigned int i=0);
952 :
953 : /**
954 : * \returns A reference to one of the system's adjoint solution
955 : * vectors, by default the one corresponding to the first qoi.
956 : */
957 : NumericVector<Number> & get_adjoint_solution(unsigned int i=0);
958 :
959 : /**
960 : * \returns A reference to one of the system's adjoint solution
961 : * vectors, by default the one corresponding to the first qoi.
962 : */
963 : const NumericVector<Number> & get_adjoint_solution(unsigned int i=0) const;
964 :
965 : /**
966 : * \returns A reference to one of the system's solution sensitivity
967 : * vectors, by default the one corresponding to the first parameter.
968 : * Creates the vector if it doesn't already exist.
969 : */
970 : NumericVector<Number> & add_sensitivity_solution(unsigned int i=0);
971 :
972 : /**
973 : * \returns A reference to one of the system's solution sensitivity
974 : * vectors, by default the one corresponding to the first parameter.
975 : */
976 : NumericVector<Number> & get_sensitivity_solution(unsigned int i=0);
977 :
978 : /**
979 : * \returns A reference to one of the system's solution sensitivity
980 : * vectors, by default the one corresponding to the first parameter.
981 : */
982 : const NumericVector<Number> & get_sensitivity_solution(unsigned int i=0) const;
983 :
984 : /**
985 : * \returns A reference to one of the system's weighted sensitivity
986 : * adjoint solution vectors, by default the one corresponding to the
987 : * first qoi.
988 : * Creates the vector if it doesn't already exist.
989 : */
990 : NumericVector<Number> & add_weighted_sensitivity_adjoint_solution(unsigned int i=0);
991 :
992 : /**
993 : * \returns A reference to one of the system's weighted sensitivity
994 : * adjoint solution vectors, by default the one corresponding to the
995 : * first qoi.
996 : */
997 : NumericVector<Number> & get_weighted_sensitivity_adjoint_solution(unsigned int i=0);
998 :
999 : /**
1000 : * \returns A reference to one of the system's weighted sensitivity
1001 : * adjoint solution vectors, by default the one corresponding to the
1002 : * first qoi.
1003 : */
1004 : const NumericVector<Number> & get_weighted_sensitivity_adjoint_solution(unsigned int i=0) const;
1005 :
1006 : /**
1007 : * \returns A reference to the solution of the last weighted
1008 : * sensitivity solve
1009 : * Creates the vector if it doesn't already exist.
1010 : */
1011 : NumericVector<Number> & add_weighted_sensitivity_solution();
1012 :
1013 : /**
1014 : * \returns A reference to the solution of the last weighted
1015 : * sensitivity solve
1016 : */
1017 : NumericVector<Number> & get_weighted_sensitivity_solution();
1018 :
1019 : /**
1020 : * \returns A reference to the solution of the last weighted
1021 : * sensitivity solve
1022 : */
1023 : const NumericVector<Number> & get_weighted_sensitivity_solution() const;
1024 :
1025 : /**
1026 : * \returns A reference to one of the system's adjoint rhs
1027 : * vectors, by default the one corresponding to the first qoi.
1028 : * Creates the vector if it doesn't already exist.
1029 : */
1030 : NumericVector<Number> & add_adjoint_rhs(unsigned int i=0);
1031 :
1032 : /**
1033 : * \returns A reference to one of the system's adjoint rhs
1034 : * vectors, by default the one corresponding to the first qoi.
1035 : * This what the user's QoI derivative code should assemble
1036 : * when setting up an adjoint problem
1037 : */
1038 : NumericVector<Number> & get_adjoint_rhs(unsigned int i=0);
1039 :
1040 : /**
1041 : * \returns A reference to one of the system's adjoint rhs
1042 : * vectors, by default the one corresponding to the first qoi.
1043 : */
1044 : const NumericVector<Number> & get_adjoint_rhs(unsigned int i=0) const;
1045 :
1046 : /**
1047 : * \returns A reference to one of the system's sensitivity rhs
1048 : * vectors, by default the one corresponding to the first parameter.
1049 : * Creates the vector if it doesn't already exist.
1050 : */
1051 : NumericVector<Number> & add_sensitivity_rhs(unsigned int i=0);
1052 :
1053 : /**
1054 : * \returns A reference to one of the system's sensitivity rhs
1055 : * vectors, by default the one corresponding to the first parameter.
1056 : * By default these vectors are built by the library, using finite
1057 : * differences, when \p assemble_residual_derivatives() is called.
1058 : *
1059 : * When assembled, this vector should hold
1060 : * -(partial R / partial p_i)
1061 : */
1062 : NumericVector<Number> & get_sensitivity_rhs(unsigned int i=0);
1063 :
1064 : /**
1065 : * \returns A reference to one of the system's sensitivity rhs
1066 : * vectors, by default the one corresponding to the first parameter.
1067 : */
1068 : const NumericVector<Number> & get_sensitivity_rhs(unsigned int i=0) const;
1069 :
1070 : /**
1071 : * \returns The number of vectors (in addition to the solution)
1072 : * handled by this system
1073 : * This is the size of the \p _vectors map
1074 : */
1075 : unsigned int n_vectors () const;
1076 :
1077 : /**
1078 : * \returns The number of matrices
1079 : * handled by this system.
1080 : * This is the size of the \p _matrices map
1081 : */
1082 : unsigned int n_matrices () const;
1083 :
1084 : /**
1085 : * \returns The number of variables in the system
1086 : */
1087 : unsigned int n_vars() const;
1088 :
1089 : /**
1090 : * \returns The number of \p VariableGroup variable groups in the system
1091 : */
1092 : unsigned int n_variable_groups() const;
1093 :
1094 : /**
1095 : * \returns The total number of scalar components in the system's
1096 : * variables. This will equal \p n_vars() in the case of all
1097 : * scalar-valued variables.
1098 : */
1099 : unsigned int n_components() const;
1100 :
1101 : /**
1102 : * \returns The number of degrees of freedom in the system
1103 : */
1104 : dof_id_type n_dofs() const;
1105 :
1106 : /**
1107 : * \returns The number of active degrees of freedom
1108 : * for this System.
1109 : */
1110 : dof_id_type n_active_dofs() const;
1111 :
1112 : /**
1113 : * \returns The total number of constrained degrees of freedom
1114 : * in the system.
1115 : */
1116 : dof_id_type n_constrained_dofs() const;
1117 :
1118 : /**
1119 : * \returns The number of constrained degrees of freedom
1120 : * on this processor.
1121 : */
1122 : dof_id_type n_local_constrained_dofs() const;
1123 :
1124 : /**
1125 : * \returns The number of degrees of freedom local
1126 : * to this processor
1127 : */
1128 : dof_id_type n_local_dofs() const;
1129 :
1130 : /**
1131 : * Adds the variable \p var to the list of variables
1132 : * for this system. If \p active_subdomains is either \p nullptr
1133 : * (the default) or points to an empty set, then it will be assumed that
1134 : * \p var has no subdomain restrictions
1135 : *
1136 : * \returns The index number for the new variable.
1137 : */
1138 : unsigned int add_variable (std::string_view var,
1139 : const FEType & type,
1140 : const std::set<subdomain_id_type> * const active_subdomains = nullptr);
1141 :
1142 : /**
1143 : * Adds the variable \p var to the list of variables
1144 : * for this system. Same as before, but assumes \p LAGRANGE
1145 : * as default value for \p FEType.family. If \p active_subdomains is either
1146 : * \p nullptr (the default) or points to an empty set, then it will be assumed
1147 : * that \p var has no subdomain restrictions
1148 : */
1149 : unsigned int add_variable (std::string_view var,
1150 : const Order order = FIRST,
1151 : const FEFamily = LAGRANGE,
1152 : const std::set<subdomain_id_type> * const active_subdomains = nullptr);
1153 :
1154 : /**
1155 : * Adds the variables \p vars to the list of variables
1156 : * for this system. If \p active_subdomains is either \p nullptr
1157 : * (the default) or points to an empty set, then it will be assumed that
1158 : * the \p vars have no subdomain restrictions
1159 : *
1160 : * \returns The index number for the last of the new variables.
1161 : */
1162 : unsigned int add_variables (const std::vector<std::string> & vars,
1163 : const FEType & type,
1164 : const std::set<subdomain_id_type> * const active_subdomains = nullptr);
1165 :
1166 : /**
1167 : * Adds variables \p vars to the list of variables
1168 : * for this system. If \p active_subdomains is either \p nullptr
1169 : * (the default) or points to an empty set, then it will be assumed that
1170 : * the \p vars have no subdomain restrictions. This API will end up
1171 : * calling \p this->add_variables(). However, we will additionally store data
1172 : * that can be leveraged by the \p DofMap to build degrees of freedom
1173 : * containers corresponding to all the variables in this variable array
1174 : *
1175 : * An 'array variable' is simply a sequence
1176 : * of contiguous variable numbers defined by pair where the first member of the pair
1177 : * is the first number in the variable sequence and the second member of the pair is
1178 : * the number of the last variable in the sequence plus one. Array variables may be
1179 : * used in tandem with variable grouping by downstream code to build optimized physics
1180 : * kernels since each variable in the array will have the same shape functions.
1181 : *
1182 : * \returns The index number for the last of the new variables.
1183 : */
1184 : unsigned int add_variable_array (const std::vector<std::string> & vars,
1185 : const FEType & type,
1186 : const std::set<subdomain_id_type> * const active_subdomains = nullptr);
1187 :
1188 : /**
1189 : * Adds the variable \p var to the list of variables
1190 : * for this system. Same as before, but assumes \p LAGRANGE
1191 : * as default value for \p FEType.family. If \p active_subdomains is either
1192 : * \p nullptr (the default) or points to an empty set, then it will be assumed that
1193 : * \p var has no subdomain restrictions
1194 : */
1195 : unsigned int add_variables (const std::vector<std::string> & vars,
1196 : const Order order = FIRST,
1197 : const FEFamily = LAGRANGE,
1198 : const std::set<subdomain_id_type> * const active_subdomains = nullptr);
1199 :
1200 : /**
1201 : * Return a constant reference to \p Variable \p var.
1202 : */
1203 : const Variable & variable (unsigned int var) const;
1204 :
1205 : /**
1206 : * Return a constant reference to \p VariableGroup \p vg.
1207 : */
1208 : const VariableGroup & variable_group (unsigned int vg) const;
1209 :
1210 : /**
1211 : * \returns \p true if a variable named \p var exists in this System
1212 : */
1213 : bool has_variable(std::string_view var) const;
1214 :
1215 : /**
1216 : * \returns The name of variable \p i.
1217 : */
1218 : const std::string & variable_name(const unsigned int i) const;
1219 :
1220 : /**
1221 : * \returns The variable number associated with
1222 : * the user-specified variable named \p var.
1223 : */
1224 : unsigned int variable_number (std::string_view var) const;
1225 :
1226 : /**
1227 : * Fills \p all_variable_numbers with all the variable numbers for the
1228 : * variables that have been added to this system.
1229 : */
1230 : void get_all_variable_numbers(std::vector<unsigned int> & all_variable_numbers) const;
1231 :
1232 : /**
1233 : * \returns An index, starting from 0 for the first component of the
1234 : * first variable, and incrementing for each component of each
1235 : * (potentially vector-valued) variable in the system in order.
1236 : * For systems with only scalar-valued variables, this will be the
1237 : * same as \p variable_number(var)
1238 : *
1239 : * Irony: currently our only non-scalar-valued variable type is
1240 : * SCALAR.
1241 : */
1242 : unsigned int variable_scalar_number (std::string_view var,
1243 : unsigned int component) const;
1244 :
1245 : /**
1246 : * \returns An index, starting from 0 for the first component of the
1247 : * first variable, and incrementing for each component of each
1248 : * (potentially vector-valued) variable in the system in order.
1249 : * For systems with only scalar-valued variables, this will be the
1250 : * same as \p var_num
1251 : *
1252 : * Irony: currently our only non-scalar-valued variable type is
1253 : * SCALAR.
1254 : */
1255 : unsigned int variable_scalar_number (unsigned int var_num,
1256 : unsigned int component) const;
1257 :
1258 :
1259 : /**
1260 : * \returns The finite element type variable number \p i.
1261 : */
1262 : const FEType & variable_type (const unsigned int i) const;
1263 :
1264 : /**
1265 : * \returns The finite element type for variable \p var.
1266 : */
1267 : const FEType & variable_type (std::string_view var) const;
1268 :
1269 : /**
1270 : * \returns \p true when \p VariableGroup structures should be
1271 : * automatically identified, \p false otherwise.
1272 : */
1273 : bool identify_variable_groups () const;
1274 :
1275 : /**
1276 : * Toggle automatic \p VariableGroup identification.
1277 : */
1278 : void identify_variable_groups (const bool);
1279 :
1280 : /**
1281 : * \returns A norm of variable \p var in the vector \p v, in the specified
1282 : * norm (e.g. L2, L_INF, H1)
1283 : */
1284 : Real calculate_norm(const NumericVector<Number> & v,
1285 : unsigned int var,
1286 : FEMNormType norm_type,
1287 : std::set<unsigned int> * skip_dimensions=nullptr) const;
1288 :
1289 : /**
1290 : * \returns A norm of the vector \p v, using \p component_norm and \p
1291 : * component_scale to choose and weight the norms of each variable.
1292 : */
1293 : Real calculate_norm(const NumericVector<Number> & v,
1294 : const SystemNorm & norm,
1295 : std::set<unsigned int> * skip_dimensions=nullptr) const;
1296 :
1297 : /**
1298 : * Reads the basic data header for this System.
1299 : */
1300 : void read_header (Xdr & io,
1301 : std::string_view version,
1302 : const bool read_header=true,
1303 : const bool read_additional_data=true,
1304 : const bool read_legacy_format=false);
1305 :
1306 : /**
1307 : * Reads additional data, namely vectors, for this System.
1308 : *
1309 : * \deprecated The ability to read XDR data files in the old (aka
1310 : * "legacy") XDR format has been deprecated for many years, this
1311 : * capability may soon disappear altogether.
1312 : */
1313 : #ifdef LIBMESH_ENABLE_DEPRECATED
1314 : void read_legacy_data (Xdr & io,
1315 : const bool read_additional_data=true);
1316 : #endif
1317 :
1318 : /**
1319 : * Reads additional data, namely vectors, for this System.
1320 : * This method may safely be called on a distributed-memory mesh.
1321 : */
1322 : template <typename ValType>
1323 : void read_serialized_data (Xdr & io,
1324 : const bool read_additional_data=true);
1325 : /**
1326 : * Non-templated version for backward compatibility.
1327 : *
1328 : * Reads additional data, namely vectors, for this System.
1329 : * This method may safely be called on a distributed-memory mesh.
1330 : */
1331 0 : void read_serialized_data (Xdr & io,
1332 : const bool read_additional_data=true)
1333 0 : { read_serialized_data<Number>(io, read_additional_data); }
1334 :
1335 : /**
1336 : * Read a number of identically distributed vectors. This method
1337 : * allows for optimization for the multiple vector case by only communicating
1338 : * the metadata once.
1339 : */
1340 : template <typename InValType>
1341 : std::size_t read_serialized_vectors (Xdr & io,
1342 : const std::vector<NumericVector<Number> *> & vectors) const;
1343 :
1344 : /**
1345 : * Non-templated version for backward compatibility.
1346 : *
1347 : * Read a number of identically distributed vectors. This method
1348 : * allows for optimization for the multiple vector case by only communicating
1349 : * the metadata once.
1350 : */
1351 8 : std::size_t read_serialized_vectors (Xdr & io,
1352 : const std::vector<NumericVector<Number> *> & vectors) const
1353 284 : { return read_serialized_vectors<Number>(io, vectors); }
1354 :
1355 : /**
1356 : * Reads additional data, namely vectors, for this System.
1357 : * This method may safely be called on a distributed-memory mesh.
1358 : * This method will read an individual file for each processor in the simulation
1359 : * where the local solution components for that processor are stored.
1360 : */
1361 : template <typename InValType>
1362 : void read_parallel_data (Xdr & io,
1363 : const bool read_additional_data);
1364 :
1365 : /**
1366 : * Non-templated version for backward compatibility.
1367 : *
1368 : * Reads additional data, namely vectors, for this System.
1369 : * This method may safely be called on a distributed-memory mesh.
1370 : * This method will read an individual file for each processor in the simulation
1371 : * where the local solution components for that processor are stored.
1372 : */
1373 : void read_parallel_data (Xdr & io,
1374 : const bool read_additional_data)
1375 : { read_parallel_data<Number>(io, read_additional_data); }
1376 :
1377 : /**
1378 : * Writes the basic data header for this System.
1379 : */
1380 : void write_header (Xdr & io,
1381 : std::string_view version,
1382 : const bool write_additional_data) const;
1383 :
1384 : /**
1385 : * Writes additional data, namely vectors, for this System.
1386 : * This method may safely be called on a distributed-memory mesh.
1387 : */
1388 : void write_serialized_data (Xdr & io,
1389 : const bool write_additional_data = true) const;
1390 :
1391 : /**
1392 : * Serialize & write a number of identically distributed vectors. This method
1393 : * allows for optimization for the multiple vector case by only communicating
1394 : * the metadata once.
1395 : */
1396 : std::size_t write_serialized_vectors (Xdr & io,
1397 : const std::vector<const NumericVector<Number> *> & vectors) const;
1398 :
1399 : /**
1400 : * Writes additional data, namely vectors, for this System.
1401 : * This method may safely be called on a distributed-memory mesh.
1402 : * This method will create an individual file for each processor in the simulation
1403 : * where the local solution components for that processor will be stored.
1404 : */
1405 : void write_parallel_data (Xdr & io,
1406 : const bool write_additional_data) const;
1407 :
1408 : /**
1409 : * \returns A string containing information about the
1410 : * system.
1411 : */
1412 : std::string get_info () const;
1413 :
1414 : /**
1415 : * Register a user function to use in initializing the system.
1416 : */
1417 : void attach_init_function (void fptr(EquationSystems & es,
1418 : const std::string & name));
1419 :
1420 : /**
1421 : * Register a user class to use to initialize the system.
1422 : *
1423 : * \note This is exclusive with the \p attach_init_function.
1424 : */
1425 : void attach_init_object (Initialization & init);
1426 :
1427 : /**
1428 : * Register a user function to use in assembling the system
1429 : * matrix and RHS.
1430 : */
1431 : void attach_assemble_function (void fptr(EquationSystems & es,
1432 : const std::string & name));
1433 :
1434 : /**
1435 : * Register a user object to use in assembling the system
1436 : * matrix and RHS.
1437 : */
1438 : void attach_assemble_object (Assembly & assemble);
1439 :
1440 : /**
1441 : * Register a user function for imposing constraints.
1442 : */
1443 : void attach_constraint_function (void fptr(EquationSystems & es,
1444 : const std::string & name));
1445 :
1446 : /**
1447 : * Register a user object for imposing constraints.
1448 : */
1449 : void attach_constraint_object (Constraint & constrain);
1450 :
1451 : /**
1452 : * \returns true if there is a user-defined constraint object
1453 : * attached to this object, false otherwise. Calling System::
1454 : * get_constraint_object() when there is no user-defined constraint
1455 : * object attached leads to either undefined behavior (dereferencing
1456 : * a nullptr) or an assert (in dbg mode) so you should call this
1457 : * function first unless you are sure there is a user-defined
1458 : * constraint object attached.
1459 : */
1460 : bool has_constraint_object () const;
1461 :
1462 : /**
1463 : * Return the user object for imposing constraints.
1464 : */
1465 : Constraint& get_constraint_object ();
1466 :
1467 : /**
1468 : * Register a user function for evaluating the quantities of interest,
1469 : * whose values should be placed in \p System::qoi
1470 : */
1471 : void attach_QOI_function (void fptr(EquationSystems & es,
1472 : const std::string & name,
1473 : const QoISet & qoi_indices));
1474 :
1475 : /**
1476 : * Register a user object for evaluating the quantities of interest,
1477 : * whose values should be placed in \p System::qoi
1478 : */
1479 : void attach_QOI_object (QOI & qoi);
1480 :
1481 : /**
1482 : * Register a user function for evaluating derivatives of a quantity
1483 : * of interest with respect to test functions, whose values should
1484 : * be placed in \p System::rhs
1485 : */
1486 : void attach_QOI_derivative (void fptr(EquationSystems & es,
1487 : const std::string & name,
1488 : const QoISet & qoi_indices,
1489 : bool include_liftfunc,
1490 : bool apply_constraints));
1491 :
1492 : /**
1493 : * Register a user object for evaluating derivatives of a quantity
1494 : * of interest with respect to test functions, whose values should
1495 : * be placed in \p System::rhs
1496 : */
1497 : void attach_QOI_derivative_object (QOIDerivative & qoi_derivative);
1498 :
1499 : /**
1500 : * Calls user's attached initialization function, or is overridden by
1501 : * the user in derived classes.
1502 : */
1503 : virtual void user_initialization ();
1504 :
1505 : /**
1506 : * Calls user's attached assembly function, or is overridden by
1507 : * the user in derived classes.
1508 : */
1509 : virtual void user_assembly ();
1510 :
1511 : /**
1512 : * Calls user's attached constraint function, or is overridden by
1513 : * the user in derived classes.
1514 : */
1515 : virtual void user_constrain ();
1516 :
1517 : /**
1518 : * Calls user's attached quantity of interest function, or is
1519 : * overridden by the user in derived classes.
1520 : */
1521 : virtual void user_QOI (const QoISet & qoi_indices);
1522 :
1523 : /**
1524 : * Calls user's attached quantity of interest derivative function,
1525 : * or is overridden by the user in derived classes.
1526 : */
1527 : virtual void user_QOI_derivative (const QoISet & qoi_indices = QoISet(),
1528 : bool include_liftfunc = true,
1529 : bool apply_constraints = true);
1530 :
1531 : /**
1532 : * Re-update the local values when the mesh has changed.
1533 : * This method takes the data updated by \p update() and
1534 : * makes it up-to-date on the current mesh.
1535 : */
1536 : virtual void re_update ();
1537 :
1538 : /**
1539 : * Restrict vectors after the mesh has coarsened
1540 : */
1541 : virtual void restrict_vectors ();
1542 :
1543 : /**
1544 : * Prolong vectors after the mesh has refined
1545 : */
1546 : virtual void prolong_vectors ();
1547 :
1548 : /// Parameters for the system. If a parameter is not provided, it should be retrieved from the EquationSystems
1549 : Parameters parameters;
1550 :
1551 : /**
1552 : * Flag which tells the system to whether or not to
1553 : * call the user assembly function during each call to solve().
1554 : * By default, every call to solve() begins with a call to the
1555 : * user assemble, so this flag is true. (For explicit systems,
1556 : * "solving" the system occurs during the assembly step, so this
1557 : * flag is always true for explicit systems.)
1558 : *
1559 : * You will only want to set this to false if you need direct
1560 : * control over when the system is assembled, and are willing to
1561 : * track the state of its assembly yourself. An example of such a
1562 : * case is an implicit system with multiple right hand sides. In
1563 : * this instance, a single assembly would likely be followed with
1564 : * multiple calls to solve.
1565 : *
1566 : * The frequency system and Newmark system have their own versions
1567 : * of this flag, called _finished_assemble, which might be able to
1568 : * be replaced with this more general concept.
1569 : */
1570 : bool assemble_before_solve;
1571 :
1572 : /**
1573 : * Avoids use of any cached data that might affect any solve result. Should
1574 : * be overridden in derived systems.
1575 : */
1576 : virtual void disable_cache ();
1577 :
1578 : /**
1579 : * A boolean to be set to true by systems using elem_fixed_solution,
1580 : * for optional use by e.g. stabilized methods.
1581 : * False by default.
1582 : *
1583 : * \note For FEMSystem users, if this variable is set to true, it
1584 : * must be before init_data() is called.
1585 : */
1586 : bool use_fixed_solution;
1587 :
1588 : /**
1589 : * A member int that can be employed to indicate increased or
1590 : * reduced quadrature order.
1591 : *
1592 : * \note For FEMSystem users, by default, when calling the
1593 : * user-defined residual functions, the FEMSystem will first set up
1594 : * an appropriate FEType::default_quadrature_rule() object for
1595 : * performing the integration. This rule will integrate elements of
1596 : * order up to 2*p+1 exactly (where p is the sum of the base FEType
1597 : * and local p refinement levels), but if additional (or reduced)
1598 : * quadrature accuracy is desired then this extra_quadrature_order
1599 : * (default 0) will be added.
1600 : */
1601 : int extra_quadrature_order;
1602 :
1603 :
1604 : //--------------------------------------------------
1605 : // The solution and solution access members
1606 :
1607 : /**
1608 : * \returns The current solution for the specified global
1609 : * DOF.
1610 : */
1611 : Number current_solution (const dof_id_type global_dof_number) const;
1612 :
1613 : /**
1614 : * Data structure to hold solution values.
1615 : */
1616 : std::unique_ptr<NumericVector<Number>> solution;
1617 :
1618 : /**
1619 : * All the values I need to compute my contribution
1620 : * to the simulation at hand. Think of this as the
1621 : * current solution with any ghost values needed from
1622 : * other processors. This vector is necessarily larger
1623 : * than the \p solution vector in the case of a parallel
1624 : * simulation. The \p update() member is used to synchronize
1625 : * the contents of the \p solution and \p current_local_solution
1626 : * vectors.
1627 : */
1628 : std::unique_ptr<NumericVector<Number>> current_local_solution;
1629 :
1630 : /**
1631 : * For time-dependent problems, this is the time t at the beginning of
1632 : * the current timestep.
1633 : *
1634 : * \note For DifferentiableSystem users: do \e not access this time
1635 : * during an assembly! Use the DiffContext::time value instead to
1636 : * get correct results.
1637 : */
1638 : Real time;
1639 :
1640 : /**
1641 : * Number of currently active quantities of interest.
1642 : */
1643 : unsigned int n_qois() const;
1644 :
1645 : #ifndef LIBMESH_ENABLE_DEPRECATED // We use accessors for these now
1646 : private:
1647 : #endif
1648 : /**
1649 : * Values of the quantities of interest. This vector needs
1650 : * to be both resized and filled by the user before any quantity of
1651 : * interest assembly is done and before any sensitivities are
1652 : * calculated.
1653 : */
1654 : std::vector<Number> qoi;
1655 :
1656 : /**
1657 : * Vector to hold error estimates for qois, either from a steady
1658 : * state calculation, or from a single unsteady solver timestep. Used
1659 : * by the library after resizing to match the size of the qoi vector.
1660 : * User code can use this for accumulating error estimates for example.
1661 : */
1662 : std::vector<Number> qoi_error_estimates;
1663 : #ifndef LIBMESH_ENABLE_DEPRECATED
1664 : public:
1665 : #endif
1666 :
1667 : /**
1668 : * Accessors for qoi and qoi_error_estimates vectors
1669 : */
1670 : void init_qois(unsigned int n_qois);
1671 :
1672 : void set_qoi(unsigned int qoi_index, Number qoi_value);
1673 : Number get_qoi_value(unsigned int qoi_index) const;
1674 :
1675 : void set_qoi(std::vector<Number> new_qoi);
1676 :
1677 : /**
1678 : * Returns a *copy* of qoi, not a reference
1679 : */
1680 : std::vector<Number> get_qoi_values() const;
1681 :
1682 : void set_qoi_error_estimate(unsigned int qoi_index, Number qoi_error_estimate);
1683 : Number get_qoi_error_estimate_value(unsigned int qoi_index) const;
1684 :
1685 : /**
1686 : * \returns The value of the solution variable \p var at the physical
1687 : * point \p p in the mesh, without knowing a priori which element
1688 : * contains \p p, using the degree of freedom coefficients in \p sol
1689 : * (or in \p current_local_solution if \p sol is left null).
1690 : *
1691 : * \note This function uses \p MeshBase::sub_point_locator(); users
1692 : * may or may not want to call \p MeshBase::clear_point_locator()
1693 : * afterward. Also, point_locator() is expensive (N log N for
1694 : * initial construction, log N for evaluations). Avoid using this
1695 : * function in any context where you are already looping over
1696 : * elements.
1697 : *
1698 : * Because the element containing \p p may lie on any processor,
1699 : * this function is parallel-only.
1700 : *
1701 : * By default this method expects the point to reside inside the domain
1702 : * and will abort if no element can be found which contains \p p. The
1703 : * optional parameter \p insist_on_success can be set to false to allow
1704 : * the method to return 0 when the point is not located.
1705 : */
1706 : Number point_value(unsigned int var, const Point & p,
1707 : const bool insist_on_success = true,
1708 : const NumericVector<Number> * sol = nullptr) const;
1709 :
1710 : /**
1711 : * \returns The value of the solution variable \p var at the physical
1712 : * point \p p contained in local Elem \p e, using the degree of
1713 : * freedom coefficients in \p sol (or in \p current_local_solution
1714 : * if \p sol is left null).
1715 : *
1716 : * This version of point_value can be run in serial, but assumes \p e is in
1717 : * the local mesh partition or is algebraically ghosted.
1718 : */
1719 : Number point_value(unsigned int var, const Point & p, const Elem & e,
1720 : const NumericVector<Number> * sol = nullptr) const;
1721 :
1722 : /**
1723 : * Calls the version of point_value() which takes a reference.
1724 : * This function exists only to prevent people from calling the
1725 : * version of point_value() that has a boolean third argument, which
1726 : * would result in unnecessary PointLocator calls.
1727 : */
1728 : Number point_value(unsigned int var, const Point & p, const Elem * e) const;
1729 :
1730 : /**
1731 : * Calls the parallel version of point_value().
1732 : * This function exists only to prevent people from accidentally
1733 : * calling the version of point_value() that has a boolean third
1734 : * argument, which would result in incorrect output.
1735 : */
1736 : Number point_value(unsigned int var, const Point & p, const NumericVector<Number> * sol) const;
1737 :
1738 : /**
1739 : * \returns The gradient of the solution variable \p var at the physical
1740 : * point \p p in the mesh, similarly to point_value.
1741 : */
1742 : Gradient point_gradient(unsigned int var, const Point & p,
1743 : const bool insist_on_success = true,
1744 : const NumericVector<Number> * sol = nullptr) const;
1745 :
1746 : /**
1747 : * \returns The gradient of the solution variable \p var at the physical
1748 : * point \p p in local Elem \p e in the mesh, similarly to point_value.
1749 : */
1750 : Gradient point_gradient(unsigned int var, const Point & p, const Elem & e,
1751 : const NumericVector<Number> * sol = nullptr) const;
1752 :
1753 : /**
1754 : * Calls the version of point_gradient() which takes a reference.
1755 : * This function exists only to prevent people from calling the
1756 : * version of point_gradient() that has a boolean third argument, which
1757 : * would result in unnecessary PointLocator calls.
1758 : */
1759 : Gradient point_gradient(unsigned int var, const Point & p, const Elem * e) const;
1760 :
1761 : /**
1762 : * Calls the parallel version of point_gradient().
1763 : * This function exists only to prevent people from accidentally
1764 : * calling the version of point_gradient() that has a boolean third
1765 : * argument, which would result in incorrect output.
1766 : */
1767 : Gradient point_gradient(unsigned int var, const Point & p, const NumericVector<Number> * sol) const;
1768 :
1769 : /**
1770 : * \returns The second derivative tensor of the solution variable \p var
1771 : * at the physical point \p p in the mesh, similarly to point_value.
1772 : */
1773 : Tensor point_hessian(unsigned int var, const Point & p,
1774 : const bool insist_on_success = true,
1775 : const NumericVector<Number> * sol = nullptr) const;
1776 :
1777 : /**
1778 : * \returns The second derivative tensor of the solution variable \p var
1779 : * at the physical point \p p in local Elem \p e in the mesh, similarly to
1780 : * point_value.
1781 : */
1782 : Tensor point_hessian(unsigned int var, const Point & p, const Elem & e,
1783 : const NumericVector<Number> * sol = nullptr) const;
1784 :
1785 : /**
1786 : * Calls the version of point_hessian() which takes a reference.
1787 : * This function exists only to prevent people from calling the
1788 : * version of point_hessian() that has a boolean third argument, which
1789 : * would result in unnecessary PointLocator calls.
1790 : */
1791 : Tensor point_hessian(unsigned int var, const Point & p, const Elem * e) const;
1792 :
1793 : /**
1794 : * Calls the parallel version of point_hessian().
1795 : * This function exists only to prevent people from accidentally
1796 : * calling the version of point_hessian() that has a boolean third
1797 : * argument, which would result in incorrect output.
1798 : */
1799 : Tensor point_hessian(unsigned int var, const Point & p, const NumericVector<Number> * sol) const;
1800 :
1801 :
1802 : /**
1803 : * Fills the std::set with the degrees of freedom on the local
1804 : * processor corresponding the the variable number passed in.
1805 : */
1806 : void local_dof_indices (const unsigned int var,
1807 : std::set<dof_id_type> & var_indices) const;
1808 :
1809 : /**
1810 : * Zeroes all dofs in \p v that correspond to variable number \p
1811 : * var_num.
1812 : */
1813 : void zero_variable (NumericVector<Number> & v, unsigned int var_num) const;
1814 :
1815 : /**
1816 : * Setter and getter functions for project_with_constraints boolean
1817 : */
1818 476 : bool get_project_with_constraints()
1819 : {
1820 16190 : return project_with_constraints;
1821 : }
1822 :
1823 952 : void set_project_with_constraints(bool _project_with_constraints)
1824 : {
1825 32856 : project_with_constraints = _project_with_constraints;
1826 952 : }
1827 :
1828 : /**
1829 : * \returns A writable reference to a boolean that determines if this system
1830 : * can be written to file or not. If set to \p true, then
1831 : * \p EquationSystems::write will ignore this system.
1832 : */
1833 11774 : bool & hide_output() { return _hide_output; }
1834 :
1835 : #ifdef LIBMESH_HAVE_METAPHYSICL
1836 : /**
1837 : * This method creates a projection matrix which corresponds to the
1838 : * operation of project_vector between old and new solution spaces.
1839 : *
1840 : * Heterogeneous Dirichlet boundary conditions are *not* taken into
1841 : * account here; if this matrix is used for prolongation (mesh
1842 : * refinement) on a side with a heterogeneous BC, the newly created
1843 : * degrees of freedom on that side will still match the coarse grid
1844 : * approximation of the BC, not the fine grid approximation.
1845 : */
1846 : void projection_matrix (SparseMatrix<Number> & proj_mat) const;
1847 : #endif // LIBMESH_HAVE_METAPHYSICL
1848 :
1849 : /**
1850 : * Adds the additional matrix \p mat_name to this system. Only
1851 : * allowed prior to \p assemble(). All additional matrices
1852 : * have the same sparsity pattern as the matrix used during
1853 : * solution. When not \p System but the user wants to
1854 : * initialize the main/system matrix, then all the additional matrices,
1855 : * if existent, have to be initialized by the user, too.
1856 : *
1857 : * This non-template method will add a derived matrix type corresponding to
1858 : * the solver package. If the user wishes to specify the matrix type to add,
1859 : * use the templated \p add_matrix method instead
1860 : *
1861 : * @param mat_name A name for the matrix
1862 : * @param type The serial/parallel/ghosted type of the matrix
1863 : * @param mat_build_type The matrix type to build
1864 : *
1865 : */
1866 : SparseMatrix<Number> & add_matrix (std::string_view mat_name,
1867 : ParallelType type = PARALLEL,
1868 : MatrixBuildType mat_build_type = MatrixBuildType::AUTOMATIC);
1869 :
1870 : /**
1871 : * Adds the additional matrix \p mat_name to this system. Only
1872 : * allowed prior to \p assemble(). All additional matrices
1873 : * have the same sparsity pattern as the matrix used during
1874 : * solution. When not \p System but the user wants to
1875 : * initialize the main/system matrix, then all the additional matrices,
1876 : * if existent, have to be initialized by the user, too.
1877 : *
1878 : * This method will create add a derived matrix of type
1879 : * \p MatrixType<Number>. One can use the non-templated \p add_matrix method to
1880 : * add a matrix corresponding to the default solver package
1881 : *
1882 : * @param mat_name A name for the matrix
1883 : * @param type The serial/parallel/ghosted type of the matrix
1884 : */
1885 : template <template <typename> class>
1886 : SparseMatrix<Number> & add_matrix (std::string_view mat_name, ParallelType = PARALLEL);
1887 :
1888 : /**
1889 : * Adds the additional matrix \p mat_name to this system. Only
1890 : * allowed prior to \p assemble(). All additional matrices
1891 : * have the same sparsity pattern as the matrix used during
1892 : * solution. When not \p System but the user wants to
1893 : * initialize the main/system matrix, then all the additional matrices,
1894 : * if existent, have to be initialized by the user, too.
1895 : *
1896 : * @param mat_name A name for the matrix
1897 : * @param matrix The matrix we are handing over the \p System for ownership
1898 : * @param type The serial/parallel/ghosted type of the matrix
1899 : *
1900 : */
1901 : SparseMatrix<Number> & add_matrix (std::string_view mat_name,
1902 : std::unique_ptr<SparseMatrix<Number>> matrix,
1903 : ParallelType type = PARALLEL);
1904 :
1905 : /**
1906 : * Removes the additional matrix \p mat_name from this system
1907 : */
1908 : void remove_matrix(std::string_view mat_name);
1909 :
1910 : /**
1911 : * \returns \p true if this \p System has a matrix associated with the
1912 : * given name, \p false otherwise.
1913 : */
1914 0 : inline bool have_matrix (std::string_view mat_name) const { return _matrices.count(mat_name); }
1915 :
1916 : /**
1917 : * \returns A const pointer to this system's additional matrix
1918 : * named \p mat_name, or \p nullptr if no matrix by that name
1919 : * exists.
1920 : */
1921 : const SparseMatrix<Number> * request_matrix (std::string_view mat_name) const;
1922 :
1923 : /**
1924 : * \returns A writable pointer to this system's additional matrix
1925 : * named \p mat_name, or \p nullptr if no matrix by that name
1926 : * exists.
1927 : */
1928 : SparseMatrix<Number> * request_matrix (std::string_view mat_name);
1929 :
1930 : /**
1931 : * \returns A const reference to this system's matrix
1932 : * named \p mat_name.
1933 : */
1934 : const SparseMatrix<Number> & get_matrix (std::string_view mat_name) const;
1935 :
1936 : /**
1937 : * \returns A writable reference to this system's matrix
1938 : * named \p mat_name.
1939 : */
1940 : SparseMatrix<Number> & get_matrix (std::string_view mat_name);
1941 :
1942 : /**
1943 : * Sets whether to use hash table matrix assembly if the matrix sub-classes support it
1944 : */
1945 : void prefer_hash_table_matrix_assembly(bool preference);
1946 :
1947 : /**
1948 : * Instructs this system to prefix solve options with its name for solvers that leverage prefixes
1949 : */
1950 0 : void prefix_with_name(bool value) { _prefix_with_name = value; }
1951 :
1952 : /**
1953 : * @returns Whether we are name prefixing
1954 : */
1955 224786 : [[nodiscard]] bool prefix_with_name() const { return _prefix_with_name; }
1956 :
1957 : /**
1958 : * @returns A prefix that may be applied to solver options. Note that this
1959 : * prefix is only used if \p prefix_with_name()
1960 : */
1961 0 : std::string prefix() const { return this->name() + "_"; }
1962 :
1963 : /**
1964 : * Request that static condensation be performed for this system
1965 : */
1966 : virtual void create_static_condensation();
1967 :
1968 : /**
1969 : * @returns Whether this system will be statically condensed
1970 : */
1971 : bool has_static_condensation() const;
1972 :
1973 : protected:
1974 :
1975 : /**
1976 : * Initializes the data for the system.
1977 : *
1978 : * \note This is called before any user-supplied initialization
1979 : * function so that all required storage will be available.
1980 : */
1981 : virtual void init_data ();
1982 :
1983 : /**
1984 : * Insertion point for adding matrices in derived classes
1985 : * before init_matrices() is called.
1986 : */
1987 220042 : virtual void add_matrices() {}
1988 :
1989 : /**
1990 : * Initializes the matrices associated with this system.
1991 : */
1992 : virtual void init_matrices ();
1993 :
1994 : /**
1995 : * \returns Whether or not matrices can still be added without
1996 : * expensive per-matrix initialization.
1997 : */
1998 272 : bool can_add_matrices() const { return !_matrices_initialized; }
1999 :
2000 : /**
2001 : * Projects the vector defined on the old mesh onto the
2002 : * new mesh.
2003 : *
2004 : * Constrain the new vector using the requested adjoint rather than
2005 : * primal constraints if is_adjoint is non-negative.
2006 : */
2007 : void project_vector (NumericVector<Number> &,
2008 : int is_adjoint = -1) const;
2009 :
2010 : /**
2011 : * Projects the vector defined on the old mesh onto the
2012 : * new mesh. The original vector is unchanged and the new vector
2013 : * is passed through the second argument.
2014 : *
2015 : * Constrain the new vector using the requested adjoint rather than
2016 : * primal constraints if is_adjoint is non-negative.
2017 : */
2018 : void project_vector (const NumericVector<Number> &,
2019 : NumericVector<Number> &,
2020 : int is_adjoint = -1) const;
2021 :
2022 : /*
2023 : * If we have e.g. a element space constrained by spline values, we
2024 : * can directly project only on the constrained basis; to get
2025 : * consistent constraining values we have to solve for them.
2026 : *
2027 : * Constrain the new vector using the requested adjoint rather than
2028 : * primal constraints if is_adjoint is non-negative.
2029 : */
2030 : void solve_for_unconstrained_dofs (NumericVector<Number> &,
2031 : int is_adjoint = -1) const;
2032 :
2033 : /**
2034 : * Whether this object should condense out constrained degrees of freedom
2035 : */
2036 140 : virtual bool condense_constrained_dofs() const { return false; }
2037 :
2038 : private:
2039 : /**
2040 : * Helper function to keep DofMap forward declarable in system.h
2041 : */
2042 : void late_matrix_init(SparseMatrix<Number> & mat,
2043 : ParallelType type);
2044 :
2045 : /**
2046 : * Finds the discrete norm for the entries in the vector
2047 : * corresponding to Dofs associated with var.
2048 : */
2049 : Real discrete_var_norm (const NumericVector<Number> & v,
2050 : unsigned int var,
2051 : FEMNormType norm_type) const;
2052 :
2053 : /**
2054 : * Reads an input vector from the stream \p io and assigns
2055 : * the values to a set of \p DofObjects. This method uses
2056 : * blocked input and is safe to call on a distributed memory-mesh.
2057 : * Unless otherwise specified, all variables are read.
2058 : *
2059 : * If an entry in \p vecs is a null pointer, the corresponding data
2060 : * is read (incrementing the file read location) but discarded.
2061 : */
2062 : template <typename iterator_type, typename InValType>
2063 : std::size_t read_serialized_blocked_dof_objects (const dof_id_type n_objects,
2064 : const iterator_type begin,
2065 : const iterator_type end,
2066 : const InValType dummy,
2067 : Xdr & io,
2068 : const std::vector<NumericVector<Number> *> & vecs,
2069 : const unsigned int var_to_read=libMesh::invalid_uint) const;
2070 :
2071 : /**
2072 : * Reads the SCALAR dofs from the stream \p io and assigns the values
2073 : * to the appropriate entries of \p vec.
2074 : *
2075 : * \returns The number of dofs read.
2076 : *
2077 : * Reads data and discards it if \p vec is a null pointer.
2078 : */
2079 : unsigned int read_SCALAR_dofs (const unsigned int var,
2080 : Xdr & io,
2081 : NumericVector<Number> * vec) const;
2082 :
2083 : /**
2084 : * Reads a vector for this System.
2085 : * This method may safely be called on a distributed-memory mesh.
2086 : *
2087 : * \returns The length of the vector read.
2088 : *
2089 : * Reads data and discards it if \p vec is a null pointer.
2090 : */
2091 : template <typename InValType>
2092 : numeric_index_type read_serialized_vector (Xdr & io,
2093 : NumericVector<Number> * vec);
2094 :
2095 : /**
2096 : * Non-templated version for backward compatibility.
2097 : *
2098 : * Reads a vector for this System.
2099 : * This method may safely be called on a distributed-memory mesh.
2100 : *
2101 : * \returns The length of the vector read.
2102 : */
2103 : numeric_index_type read_serialized_vector (Xdr & io,
2104 : NumericVector<Number> & vec)
2105 : { return read_serialized_vector<Number>(io, &vec); }
2106 :
2107 : /**
2108 : * Writes an output vector to the stream \p io for a set of \p DofObjects.
2109 : * This method uses blocked output and is safe to call on a distributed memory-mesh.
2110 : *
2111 : * \returns The number of values written
2112 : */
2113 : template <typename iterator_type>
2114 : std::size_t write_serialized_blocked_dof_objects (const std::vector<const NumericVector<Number> *> & vecs,
2115 : const dof_id_type n_objects,
2116 : const iterator_type begin,
2117 : const iterator_type end,
2118 : Xdr & io,
2119 : const unsigned int var_to_write=libMesh::invalid_uint) const;
2120 :
2121 : /**
2122 : * Writes the SCALAR dofs associated with var to the stream \p io.
2123 : *
2124 : * \returns The number of values written.
2125 : */
2126 : unsigned int write_SCALAR_dofs (const NumericVector<Number> & vec,
2127 : const unsigned int var,
2128 : Xdr & io) const;
2129 :
2130 : /**
2131 : * Writes a vector for this System.
2132 : * This method may safely be called on a distributed-memory mesh.
2133 : *
2134 : * \returns The number of values written.
2135 : */
2136 : dof_id_type write_serialized_vector (Xdr & io,
2137 : const NumericVector<Number> & vec) const;
2138 :
2139 : /**
2140 : * Function that initializes the system.
2141 : */
2142 : void (* _init_system_function) (EquationSystems & es,
2143 : const std::string & name);
2144 :
2145 : /**
2146 : * Object that initializes the system.
2147 : */
2148 : Initialization * _init_system_object;
2149 :
2150 : /**
2151 : * Function that assembles the system.
2152 : */
2153 : void (* _assemble_system_function) (EquationSystems & es,
2154 : const std::string & name);
2155 :
2156 : /**
2157 : * Object that assembles the system.
2158 : */
2159 : Assembly * _assemble_system_object;
2160 :
2161 : /**
2162 : * Function to impose constraints.
2163 : */
2164 : void (* _constrain_system_function) (EquationSystems & es,
2165 : const std::string & name);
2166 :
2167 : /**
2168 : * Object that constrains the system.
2169 : */
2170 : Constraint * _constrain_system_object;
2171 :
2172 : /**
2173 : * Function to evaluate quantity of interest
2174 : */
2175 : void (* _qoi_evaluate_function) (EquationSystems & es,
2176 : const std::string & name,
2177 : const QoISet & qoi_indices);
2178 :
2179 : /**
2180 : * Object to compute quantities of interest.
2181 : */
2182 : QOI * _qoi_evaluate_object;
2183 :
2184 : /**
2185 : * Function to evaluate quantity of interest derivative
2186 : */
2187 : void (* _qoi_evaluate_derivative_function) (EquationSystems & es,
2188 : const std::string & name,
2189 : const QoISet & qoi_indices,
2190 : bool include_liftfunc,
2191 : bool apply_constraints);
2192 :
2193 : /**
2194 : * Object to compute derivatives of quantities of interest.
2195 : */
2196 : QOIDerivative * _qoi_evaluate_derivative_object;
2197 :
2198 : /**
2199 : * Data structure describing the relationship between
2200 : * nodes, variables, etc... and degrees of freedom.
2201 : */
2202 : std::unique_ptr<DofMap> _dof_map;
2203 :
2204 : /**
2205 : * Constant reference to the \p EquationSystems object
2206 : * used for the simulation.
2207 : */
2208 : EquationSystems & _equation_systems;
2209 :
2210 : /**
2211 : * Constant reference to the \p mesh data structure used
2212 : * for the simulation.
2213 : */
2214 : MeshBase & _mesh;
2215 :
2216 : /**
2217 : * A name associated with this system.
2218 : */
2219 : const std::string _sys_name;
2220 :
2221 : /**
2222 : * The number associated with this system
2223 : */
2224 : const unsigned int _sys_number;
2225 :
2226 : /**
2227 : * Flag stating if the system is active or not.
2228 : */
2229 : bool _active;
2230 :
2231 : /**
2232 : * Some systems need an arbitrary number of vectors.
2233 : * This map allows names to be associated with arbitrary
2234 : * vectors. All the vectors in this map will be distributed
2235 : * in the same way as the solution vector.
2236 : */
2237 : std::map<std::string, std::unique_ptr<NumericVector<Number>>, std::less<>> _vectors;
2238 :
2239 : /**
2240 : * Holds true if a vector by that name should be projected
2241 : * onto a changed grid, false if it should be zeroed.
2242 : */
2243 : std::map<std::string, bool, std::less<>> _vector_projections;
2244 :
2245 : /**
2246 : * Holds non-negative if a vector by that name should be projected
2247 : * using adjoint constraints/BCs, -1 if primal
2248 : */
2249 : std::map<std::string, int, std::less<>> _vector_is_adjoint;
2250 :
2251 : /**
2252 : * Some systems need an arbitrary number of matrices.
2253 : */
2254 : std::map<std::string, std::unique_ptr<SparseMatrix<Number>>, std::less<>> _matrices;
2255 :
2256 : /**
2257 : * Holds the types of the matrices
2258 : */
2259 : std::map<std::string, ParallelType, std::less<>> _matrix_types;
2260 :
2261 : /**
2262 : * \p false when additional matrices being added require initialization, \p true otherwise.
2263 : */
2264 : bool _matrices_initialized;
2265 :
2266 : /**
2267 : * Holds true if the solution vector should be projected
2268 : * onto a changed grid, false if it should be zeroed.
2269 : * This is true by default.
2270 : */
2271 : bool _solution_projection;
2272 :
2273 : /**
2274 : * Holds true if the components of more advanced system types (e.g.
2275 : * system matrices) should not be initialized.
2276 : */
2277 : bool _basic_system_only;
2278 :
2279 : /**
2280 : * \p true when additional vectors and variables do not require
2281 : * immediate initialization, \p false otherwise.
2282 : */
2283 : bool _is_initialized;
2284 :
2285 : /**
2286 : * This flag is used only when *reading* in a system from file.
2287 : * Based on the system header, it keeps track of how many
2288 : * additional vectors were actually written for this file.
2289 : */
2290 : unsigned int _additional_data_written;
2291 :
2292 : /**
2293 : * This vector is used only when *reading* in a system from file.
2294 : * Based on the system header, it keeps track of any index remapping
2295 : * between variable names in the data file and variable names in the
2296 : * already-constructed system. I.e. if we have a system with
2297 : * variables "A1", "A2", "B1", and "B2", but we read in a data file with
2298 : * only "A1" and "B1" defined, then we don't want to try and read in
2299 : * A2 or B2, and we don't want to assign A1 and B1 values to
2300 : * different dof indices.
2301 : */
2302 : std::vector<unsigned int> _written_var_indices;
2303 :
2304 : /**
2305 : * Has the adjoint problem already been solved? If the user sets
2306 : * \p adjoint_already_solved to \p true, we won't waste time solving
2307 : * it again.
2308 : */
2309 : bool adjoint_already_solved;
2310 :
2311 : /**
2312 : * Are we allowed to write this system to file? If \p _hide_output is
2313 : * \p true, then \p EquationSystems::write will ignore this system.
2314 : */
2315 : bool _hide_output;
2316 :
2317 : /**
2318 : * Do we want to apply constraints while projecting vectors ?
2319 : */
2320 : bool project_with_constraints;
2321 :
2322 : /**
2323 : * Whether to use hash table matrix assembly if the matrix sub-classes support it
2324 : */
2325 : bool _prefer_hash_table_matrix_assembly;
2326 :
2327 : /**
2328 : * Whether any of our matrices require an initial sparsity pattern computation in order to determine preallocation
2329 : */
2330 : bool _require_sparsity_pattern;
2331 :
2332 : /**
2333 : * Whether we are name prefixing solver options
2334 : */
2335 : bool _prefix_with_name;
2336 : };
2337 :
2338 :
2339 :
2340 : // ------------------------------------------------------------
2341 : // System inline methods
2342 : inline
2343 681881 : const std::string & System::name() const
2344 : {
2345 14330898 : return _sys_name;
2346 : }
2347 :
2348 :
2349 :
2350 : inline
2351 651289 : unsigned int System::number() const
2352 : {
2353 12252659 : return _sys_number;
2354 : }
2355 :
2356 :
2357 :
2358 : inline
2359 837121 : const MeshBase & System::get_mesh() const
2360 : {
2361 17003306 : return _mesh;
2362 : }
2363 :
2364 :
2365 :
2366 : inline
2367 28734 : MeshBase & System::get_mesh()
2368 : {
2369 1740402 : return _mesh;
2370 : }
2371 :
2372 :
2373 :
2374 : inline
2375 76024217 : const DofMap & System::get_dof_map() const
2376 : {
2377 76024217 : return *_dof_map;
2378 : }
2379 :
2380 :
2381 :
2382 : inline
2383 326446 : DofMap & System::get_dof_map()
2384 : {
2385 326446 : return *_dof_map;
2386 : }
2387 :
2388 :
2389 :
2390 : inline
2391 : bool System::active() const
2392 : {
2393 : return _active;
2394 : }
2395 :
2396 :
2397 :
2398 : inline
2399 : void System::activate ()
2400 : {
2401 : _active = true;
2402 : }
2403 :
2404 :
2405 :
2406 : inline
2407 : void System::deactivate ()
2408 : {
2409 : _active = false;
2410 : }
2411 :
2412 :
2413 :
2414 : inline
2415 22324 : bool System::is_initialized () const
2416 : {
2417 44189 : return _is_initialized;
2418 : }
2419 :
2420 :
2421 :
2422 : inline
2423 0 : void System::set_basic_system_only ()
2424 : {
2425 0 : _basic_system_only = true;
2426 0 : }
2427 :
2428 :
2429 :
2430 : inline
2431 : unsigned int
2432 0 : System::variable_scalar_number (std::string_view var,
2433 : unsigned int component) const
2434 : {
2435 0 : return variable_scalar_number(this->variable_number(var), component);
2436 : }
2437 :
2438 :
2439 :
2440 : inline
2441 22218 : dof_id_type System::n_active_dofs() const
2442 : {
2443 22948 : return this->n_dofs() - this->n_constrained_dofs();
2444 : }
2445 :
2446 :
2447 :
2448 : inline
2449 : bool System::have_vector (std::string_view vec_name) const
2450 : {
2451 : return (_vectors.count(vec_name));
2452 : }
2453 :
2454 :
2455 :
2456 : inline
2457 1666 : unsigned int System::n_vectors () const
2458 : {
2459 1666 : return cast_int<unsigned int>(_vectors.size());
2460 : }
2461 :
2462 : inline
2463 1440 : System::vectors_iterator System::vectors_begin ()
2464 : {
2465 1440 : return _vectors.begin();
2466 : }
2467 :
2468 : inline
2469 1892 : System::const_vectors_iterator System::vectors_begin () const
2470 : {
2471 1892 : return _vectors.begin();
2472 : }
2473 :
2474 : inline
2475 1440 : System::vectors_iterator System::vectors_end ()
2476 : {
2477 1440 : return _vectors.end();
2478 : }
2479 :
2480 : inline
2481 3784 : System::const_vectors_iterator System::vectors_end () const
2482 : {
2483 3784 : return _vectors.end();
2484 : }
2485 :
2486 : inline
2487 : System::matrices_iterator System::matrices_begin ()
2488 : {
2489 : return _matrices.begin();
2490 : }
2491 :
2492 : inline
2493 : System::const_matrices_iterator System::matrices_begin () const
2494 : {
2495 : return _matrices.begin();
2496 : }
2497 :
2498 : inline
2499 : System::matrices_iterator System::matrices_end ()
2500 : {
2501 : return _matrices.end();
2502 : }
2503 :
2504 : inline
2505 : System::const_matrices_iterator System::matrices_end () const
2506 : {
2507 : return _matrices.end();
2508 : }
2509 :
2510 : inline
2511 0 : void System::assemble_residual_derivatives (const ParameterVector &)
2512 : {
2513 0 : libmesh_not_implemented();
2514 : }
2515 :
2516 : inline
2517 0 : void System::disable_cache () { assemble_before_solve = true; }
2518 :
2519 : inline
2520 1669214 : unsigned int System::n_qois() const
2521 : {
2522 : #ifndef LIBMESH_ENABLE_DEPRECATED
2523 : libmesh_assert_equal_to(this->qoi.size(), this->qoi_error_estimates.size());
2524 : #endif
2525 :
2526 8130634 : return cast_int<unsigned int>(this->qoi.size());
2527 : }
2528 :
2529 : inline
2530 : std::pair<unsigned int, Real>
2531 0 : System::sensitivity_solve (const ParameterVector &)
2532 : {
2533 0 : libmesh_not_implemented();
2534 : }
2535 :
2536 : inline
2537 : std::pair<unsigned int, Real>
2538 0 : System::weighted_sensitivity_solve (const ParameterVector &,
2539 : const ParameterVector &)
2540 : {
2541 0 : libmesh_not_implemented();
2542 : }
2543 :
2544 : inline
2545 : std::pair<unsigned int, Real>
2546 0 : System::adjoint_solve (const QoISet &)
2547 : {
2548 0 : libmesh_not_implemented();
2549 : }
2550 :
2551 : inline
2552 : std::pair<unsigned int, Real>
2553 0 : System::weighted_sensitivity_adjoint_solve (const ParameterVector &,
2554 : const ParameterVector &,
2555 : const QoISet &)
2556 : {
2557 0 : libmesh_not_implemented();
2558 : }
2559 :
2560 : inline
2561 : void
2562 0 : System::adjoint_qoi_parameter_sensitivity (const QoISet &,
2563 : const ParameterVector &,
2564 : SensitivityData &)
2565 : {
2566 0 : libmesh_not_implemented();
2567 : }
2568 :
2569 : inline
2570 : void
2571 0 : System::forward_qoi_parameter_sensitivity (const QoISet &,
2572 : const ParameterVector &,
2573 : SensitivityData &)
2574 : {
2575 0 : libmesh_not_implemented();
2576 : }
2577 :
2578 : inline
2579 : void
2580 0 : System::qoi_parameter_hessian(const QoISet &,
2581 : const ParameterVector &,
2582 : SensitivityData &)
2583 : {
2584 0 : libmesh_not_implemented();
2585 : }
2586 :
2587 : inline
2588 : void
2589 0 : System::qoi_parameter_hessian_vector_product(const QoISet &,
2590 : const ParameterVector &,
2591 : const ParameterVector &,
2592 : SensitivityData &)
2593 : {
2594 0 : libmesh_not_implemented();
2595 : }
2596 :
2597 : inline
2598 1058 : unsigned int System::n_matrices () const
2599 : {
2600 1058 : return cast_int<unsigned int>(_matrices.size());
2601 : }
2602 :
2603 : template <template <typename> class MatrixType>
2604 : inline
2605 : SparseMatrix<Number> &
2606 : System::add_matrix (std::string_view mat_name,
2607 : const ParallelType type)
2608 : {
2609 : // Return the matrix if it is already there.
2610 : auto it = this->_matrices.find(mat_name);
2611 : if (it != this->_matrices.end())
2612 : return *it->second;
2613 :
2614 : // Otherwise build the matrix to return.
2615 : auto pr = _matrices.emplace(mat_name, std::make_unique<MatrixType<Number>>(this->comm()));
2616 : _matrix_types.emplace(mat_name, type);
2617 :
2618 : SparseMatrix<Number> & mat = *(pr.first->second);
2619 :
2620 : // Initialize it first if we've already initialized the others.
2621 : this->late_matrix_init(mat, type);
2622 :
2623 : return mat;
2624 : }
2625 :
2626 : inline void
2627 : System::prefer_hash_table_matrix_assembly(const bool preference)
2628 : {
2629 : libmesh_error_msg_if(
2630 : _matrices_initialized,
2631 : "System::prefer_hash_table_matrix_assembly() should be called before matrices are initialized");
2632 : _prefer_hash_table_matrix_assembly = preference;
2633 : }
2634 : } // namespace libMesh
2635 :
2636 : #endif // LIBMESH_SYSTEM_H
|
