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