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