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