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 25384 : { 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 12 : 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 589011 : const EquationSystems & get_equation_systems() const { return _equation_systems; }
768 :
769 : /**
770 : * \returns A reference to this system's parent EquationSystems object.
771 : */
772 223413 : 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 284 : { 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 16190 : return project_with_constraints;
1840 : }
1841 :
1842 952 : void set_project_with_constraints(bool _project_with_constraints)
1843 : {
1844 32856 : 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 12990 : 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 231706 : [[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 0 : 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 : * If we have e.g. a element space constrained by spline values, we
1994 : * can directly project functions only on the constrained basis; to
1995 : * get consistent constraining values we have to solve for them.
1996 : *
1997 : * Constrain the new vector using the requested adjoint rather than
1998 : * primal constraints if is_adjoint is non-negative.
1999 : */
2000 : void solve_for_unconstrained_dofs (NumericVector<Number> &,
2001 : int is_adjoint = -1) const;
2002 :
2003 : protected:
2004 :
2005 : /**
2006 : * Initializes the data for the system.
2007 : *
2008 : * \note This is called before any user-supplied initialization
2009 : * function so that all required storage will be available.
2010 : */
2011 : virtual void init_data ();
2012 :
2013 : /**
2014 : * Insertion point for adding matrices in derived classes
2015 : * before init_matrices() is called.
2016 : */
2017 223124 : virtual void add_matrices() {}
2018 :
2019 : /**
2020 : * Initializes the matrices associated with this system.
2021 : */
2022 : virtual void init_matrices ();
2023 :
2024 : /**
2025 : * \returns Whether or not matrices can still be added without
2026 : * expensive per-matrix initialization.
2027 : */
2028 272 : bool can_add_matrices() const { return !_matrices_initialized; }
2029 :
2030 : /**
2031 : * Projects the vector defined on the old mesh onto the
2032 : * new mesh.
2033 : *
2034 : * Constrain the new vector using the requested adjoint rather than
2035 : * primal constraints if is_adjoint is non-negative.
2036 : */
2037 : void project_vector (NumericVector<Number> &,
2038 : int is_adjoint = -1,
2039 : std::optional<ConstElemRange> active_local_range = std::nullopt,
2040 : std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
2041 :
2042 : /**
2043 : * Projects the vector defined on the old mesh onto the
2044 : * new mesh. The original vector is unchanged and the new vector
2045 : * is passed through the second argument.
2046 : *
2047 : * Constrain the new vector using the requested adjoint rather than
2048 : * primal constraints if is_adjoint is non-negative.
2049 : */
2050 : void project_vector (const NumericVector<Number> &,
2051 : NumericVector<Number> &,
2052 : int is_adjoint = -1,
2053 : std::optional<ConstElemRange> active_local_range = std::nullopt,
2054 : std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
2055 :
2056 : /**
2057 : * Whether this object should condense out constrained degrees of freedom
2058 : */
2059 140 : virtual bool condense_constrained_dofs() const { return false; }
2060 :
2061 : private:
2062 : /**
2063 : * Helper function to keep DofMap forward declarable in system.h
2064 : */
2065 : void late_matrix_init(SparseMatrix<Number> & mat,
2066 : ParallelType type);
2067 :
2068 : /**
2069 : * Finds the discrete norm for the entries in the vector
2070 : * corresponding to Dofs associated with var.
2071 : */
2072 : Real discrete_var_norm (const NumericVector<Number> & v,
2073 : unsigned int var,
2074 : FEMNormType norm_type) const;
2075 :
2076 : /**
2077 : * Reads an input vector from the stream \p io and assigns
2078 : * the values to a set of \p DofObjects. This method uses
2079 : * blocked input and is safe to call on a distributed memory-mesh.
2080 : * Unless otherwise specified, all variables are read.
2081 : *
2082 : * If an entry in \p vecs is a null pointer, the corresponding data
2083 : * is read (incrementing the file read location) but discarded.
2084 : */
2085 : template <typename iterator_type, typename InValType>
2086 : std::size_t read_serialized_blocked_dof_objects (const dof_id_type n_objects,
2087 : const iterator_type begin,
2088 : const iterator_type end,
2089 : const InValType dummy,
2090 : Xdr & io,
2091 : const std::vector<NumericVector<Number> *> & vecs,
2092 : const unsigned int var_to_read=libMesh::invalid_uint) const;
2093 :
2094 : /**
2095 : * Reads the SCALAR dofs from the stream \p io and assigns the values
2096 : * to the appropriate entries of \p vec.
2097 : *
2098 : * \returns The number of dofs read.
2099 : *
2100 : * Reads data and discards it if \p vec is a null pointer.
2101 : */
2102 : unsigned int read_SCALAR_dofs (const unsigned int var,
2103 : Xdr & io,
2104 : NumericVector<Number> * vec) const;
2105 :
2106 : /**
2107 : * Reads a vector for this System.
2108 : * This method may safely be called on a distributed-memory mesh.
2109 : *
2110 : * \returns The length of the vector read.
2111 : *
2112 : * Reads data and discards it if \p vec is a null pointer.
2113 : */
2114 : template <typename InValType>
2115 : numeric_index_type read_serialized_vector (Xdr & io,
2116 : NumericVector<Number> * vec);
2117 :
2118 : /**
2119 : * Non-templated version for backward compatibility.
2120 : *
2121 : * Reads a vector for this System.
2122 : * This method may safely be called on a distributed-memory mesh.
2123 : *
2124 : * \returns The length of the vector read.
2125 : */
2126 : numeric_index_type read_serialized_vector (Xdr & io,
2127 : NumericVector<Number> & vec)
2128 : { return read_serialized_vector<Number>(io, &vec); }
2129 :
2130 : /**
2131 : * Writes an output vector to the stream \p io for a set of \p DofObjects.
2132 : * This method uses blocked output and is safe to call on a distributed memory-mesh.
2133 : *
2134 : * \returns The number of values written
2135 : */
2136 : template <typename iterator_type>
2137 : std::size_t write_serialized_blocked_dof_objects (const std::vector<const NumericVector<Number> *> & vecs,
2138 : const dof_id_type n_objects,
2139 : const iterator_type begin,
2140 : const iterator_type end,
2141 : Xdr & io,
2142 : const unsigned int var_to_write=libMesh::invalid_uint) const;
2143 :
2144 : /**
2145 : * Writes the SCALAR dofs associated with var to the stream \p io.
2146 : *
2147 : * \returns The number of values written.
2148 : */
2149 : unsigned int write_SCALAR_dofs (const NumericVector<Number> & vec,
2150 : const unsigned int var,
2151 : Xdr & io) const;
2152 :
2153 : /**
2154 : * Writes a vector for this System.
2155 : * This method may safely be called on a distributed-memory mesh.
2156 : *
2157 : * \returns The number of values written.
2158 : */
2159 : dof_id_type write_serialized_vector (Xdr & io,
2160 : const NumericVector<Number> & vec) const;
2161 :
2162 : /**
2163 : * Function that initializes the system.
2164 : */
2165 : void (* _init_system_function) (EquationSystems & es,
2166 : const std::string & name);
2167 :
2168 : /**
2169 : * Object that initializes the system.
2170 : */
2171 : Initialization * _init_system_object;
2172 :
2173 : /**
2174 : * Function that assembles the system.
2175 : */
2176 : void (* _assemble_system_function) (EquationSystems & es,
2177 : const std::string & name);
2178 :
2179 : /**
2180 : * Object that assembles the system.
2181 : */
2182 : Assembly * _assemble_system_object;
2183 :
2184 : /**
2185 : * Function to impose constraints.
2186 : */
2187 : void (* _constrain_system_function) (EquationSystems & es,
2188 : const std::string & name);
2189 :
2190 : /**
2191 : * Object that constrains the system.
2192 : */
2193 : Constraint * _constrain_system_object;
2194 :
2195 : /**
2196 : * Function to evaluate quantity of interest
2197 : */
2198 : void (* _qoi_evaluate_function) (EquationSystems & es,
2199 : const std::string & name,
2200 : const QoISet & qoi_indices);
2201 :
2202 : /**
2203 : * Object to compute quantities of interest.
2204 : */
2205 : QOI * _qoi_evaluate_object;
2206 :
2207 : /**
2208 : * Function to evaluate quantity of interest derivative
2209 : */
2210 : void (* _qoi_evaluate_derivative_function) (EquationSystems & es,
2211 : const std::string & name,
2212 : const QoISet & qoi_indices,
2213 : bool include_liftfunc,
2214 : bool apply_constraints);
2215 :
2216 : /**
2217 : * Object to compute derivatives of quantities of interest.
2218 : */
2219 : QOIDerivative * _qoi_evaluate_derivative_object;
2220 :
2221 : /**
2222 : * Data structure describing the relationship between
2223 : * nodes, variables, etc... and degrees of freedom.
2224 : */
2225 : std::unique_ptr<DofMap> _dof_map;
2226 :
2227 : /**
2228 : * Constant reference to the \p EquationSystems object
2229 : * used for the simulation.
2230 : */
2231 : EquationSystems & _equation_systems;
2232 :
2233 : /**
2234 : * Constant reference to the \p mesh data structure used
2235 : * for the simulation.
2236 : */
2237 : MeshBase & _mesh;
2238 :
2239 : /**
2240 : * A name associated with this system.
2241 : */
2242 : const std::string _sys_name;
2243 :
2244 : /**
2245 : * The number associated with this system
2246 : */
2247 : const unsigned int _sys_number;
2248 :
2249 : /**
2250 : * Flag stating if the system is active or not.
2251 : */
2252 : bool _active;
2253 :
2254 : /**
2255 : * Some systems need an arbitrary number of vectors.
2256 : * This map allows names to be associated with arbitrary
2257 : * vectors. All the vectors in this map will be distributed
2258 : * in the same way as the solution vector.
2259 : */
2260 : std::map<std::string, std::unique_ptr<NumericVector<Number>>, std::less<>> _vectors;
2261 :
2262 : /**
2263 : * Holds true if a vector by that name should be projected
2264 : * onto a changed grid, false if it should be zeroed.
2265 : */
2266 : std::map<std::string, bool, std::less<>> _vector_projections;
2267 :
2268 : /**
2269 : * Holds non-negative if a vector by that name should be projected
2270 : * using adjoint constraints/BCs, -1 if primal
2271 : */
2272 : std::map<std::string, int, std::less<>> _vector_is_adjoint;
2273 :
2274 : /**
2275 : * Some systems need an arbitrary number of matrices.
2276 : */
2277 : std::map<std::string, std::unique_ptr<SparseMatrix<Number>>, std::less<>> _matrices;
2278 :
2279 : /**
2280 : * Holds the types of the matrices
2281 : */
2282 : std::map<std::string, ParallelType, std::less<>> _matrix_types;
2283 :
2284 : /**
2285 : * \p false when additional matrices being added require initialization, \p true otherwise.
2286 : */
2287 : bool _matrices_initialized;
2288 :
2289 : /**
2290 : * Holds true if the solution vector should be projected
2291 : * onto a changed grid, false if it should be zeroed.
2292 : * This is true by default.
2293 : */
2294 : bool _solution_projection;
2295 :
2296 : /**
2297 : * Holds true if the components of more advanced system types (e.g.
2298 : * system matrices) should not be initialized.
2299 : */
2300 : bool _basic_system_only;
2301 :
2302 : /**
2303 : * \p true when additional vectors and variables do not require
2304 : * immediate initialization, \p false otherwise.
2305 : */
2306 : bool _is_initialized;
2307 :
2308 : /**
2309 : * This flag is used only when *reading* in a system from file.
2310 : * Based on the system header, it keeps track of how many
2311 : * additional vectors were actually written for this file.
2312 : */
2313 : unsigned int _additional_data_written;
2314 :
2315 : /**
2316 : * This vector is used only when *reading* in a system from file.
2317 : * Based on the system header, it keeps track of any index remapping
2318 : * between variable names in the data file and variable names in the
2319 : * already-constructed system. I.e. if we have a system with
2320 : * variables "A1", "A2", "B1", and "B2", but we read in a data file with
2321 : * only "A1" and "B1" defined, then we don't want to try and read in
2322 : * A2 or B2, and we don't want to assign A1 and B1 values to
2323 : * different dof indices.
2324 : */
2325 : std::vector<unsigned int> _written_var_indices;
2326 :
2327 : /**
2328 : * Has the adjoint problem already been solved? If the user sets
2329 : * \p adjoint_already_solved to \p true, we won't waste time solving
2330 : * it again.
2331 : */
2332 : bool adjoint_already_solved;
2333 :
2334 : /**
2335 : * Are we allowed to write this system to file? If \p _hide_output is
2336 : * \p true, then \p EquationSystems::write will ignore this system.
2337 : */
2338 : bool _hide_output;
2339 :
2340 : /**
2341 : * Do we want to apply constraints while projecting vectors ?
2342 : */
2343 : bool project_with_constraints;
2344 :
2345 : /**
2346 : * Whether to use hash table matrix assembly if the matrix sub-classes support it
2347 : */
2348 : bool _prefer_hash_table_matrix_assembly;
2349 :
2350 : /**
2351 : * Whether any of our matrices require an initial sparsity pattern computation in order to determine preallocation
2352 : */
2353 : bool _require_sparsity_pattern;
2354 :
2355 : /**
2356 : * Whether we are name prefixing solver options
2357 : */
2358 : bool _prefix_with_name;
2359 :
2360 : /**
2361 : * Values of the quantities of interest. This vector needs to be
2362 : * both resized and filled by the user before any quantity of
2363 : * interest assembly is done and before any sensitivities are
2364 : * calculated. Use the get_qoi_values() accessor to get these
2365 : * values.
2366 : */
2367 : std::vector<Number> _qoi;
2368 :
2369 : /**
2370 : * Vector to hold error estimates for qois, either from a steady
2371 : * state calculation, or from a single unsteady solver timestep. Used
2372 : * by the library after resizing to match the size of the qoi vector.
2373 : * User code can use this for accumulating error estimates for example.
2374 : * Use the set_qoi_error_estimate()/get_qoi_error_estimate_value()
2375 : * accessors to set/get these values.
2376 : */
2377 : std::vector<Number> _qoi_error_estimates;
2378 : };
2379 :
2380 :
2381 :
2382 : // ------------------------------------------------------------
2383 : // System inline methods
2384 : inline
2385 695075 : const std::string & System::name() const
2386 : {
2387 14486669 : return _sys_name;
2388 : }
2389 :
2390 :
2391 :
2392 : inline
2393 702397 : unsigned int System::number() const
2394 : {
2395 10945268 : return _sys_number;
2396 : }
2397 :
2398 :
2399 :
2400 : inline
2401 868499 : const MeshBase & System::get_mesh() const
2402 : {
2403 17359227 : return _mesh;
2404 : }
2405 :
2406 :
2407 :
2408 : inline
2409 28670 : MeshBase & System::get_mesh()
2410 : {
2411 1375572 : return _mesh;
2412 : }
2413 :
2414 :
2415 :
2416 : inline
2417 87561349 : const DofMap & System::get_dof_map() const
2418 : {
2419 87561349 : return *_dof_map;
2420 : }
2421 :
2422 :
2423 :
2424 : inline
2425 314142 : DofMap & System::get_dof_map()
2426 : {
2427 314142 : return *_dof_map;
2428 : }
2429 :
2430 :
2431 :
2432 : inline
2433 : bool System::active() const
2434 : {
2435 : return _active;
2436 : }
2437 :
2438 :
2439 :
2440 : inline
2441 : void System::activate ()
2442 : {
2443 : _active = true;
2444 : }
2445 :
2446 :
2447 :
2448 : inline
2449 : void System::deactivate ()
2450 : {
2451 : _active = false;
2452 : }
2453 :
2454 :
2455 :
2456 : inline
2457 23376 : bool System::is_initialized () const
2458 : {
2459 72235 : return _is_initialized;
2460 : }
2461 :
2462 :
2463 :
2464 : inline
2465 0 : void System::set_basic_system_only ()
2466 : {
2467 0 : _basic_system_only = true;
2468 0 : }
2469 :
2470 :
2471 :
2472 : inline
2473 : unsigned int
2474 0 : System::variable_scalar_number (std::string_view var,
2475 : unsigned int component) const
2476 : {
2477 0 : return variable_scalar_number(this->variable_number(var), component);
2478 : }
2479 :
2480 :
2481 :
2482 : inline
2483 22218 : dof_id_type System::n_active_dofs() const
2484 : {
2485 22948 : return this->n_dofs() - this->n_constrained_dofs();
2486 : }
2487 :
2488 :
2489 :
2490 : inline
2491 : bool System::have_vector (std::string_view vec_name) const
2492 : {
2493 : return (_vectors.count(vec_name));
2494 : }
2495 :
2496 :
2497 :
2498 : inline
2499 1664 : unsigned int System::n_vectors () const
2500 : {
2501 1664 : return cast_int<unsigned int>(_vectors.size());
2502 : }
2503 :
2504 : inline
2505 1440 : System::vectors_iterator System::vectors_begin ()
2506 : {
2507 1440 : return _vectors.begin();
2508 : }
2509 :
2510 : inline
2511 1892 : System::const_vectors_iterator System::vectors_begin () const
2512 : {
2513 1892 : return _vectors.begin();
2514 : }
2515 :
2516 : inline
2517 1440 : System::vectors_iterator System::vectors_end ()
2518 : {
2519 1440 : return _vectors.end();
2520 : }
2521 :
2522 : inline
2523 3784 : System::const_vectors_iterator System::vectors_end () const
2524 : {
2525 3784 : return _vectors.end();
2526 : }
2527 :
2528 : inline
2529 : System::matrices_iterator System::matrices_begin ()
2530 : {
2531 : return _matrices.begin();
2532 : }
2533 :
2534 : inline
2535 : System::const_matrices_iterator System::matrices_begin () const
2536 : {
2537 : return _matrices.begin();
2538 : }
2539 :
2540 : inline
2541 : System::matrices_iterator System::matrices_end ()
2542 : {
2543 : return _matrices.end();
2544 : }
2545 :
2546 : inline
2547 : System::const_matrices_iterator System::matrices_end () const
2548 : {
2549 : return _matrices.end();
2550 : }
2551 :
2552 : inline
2553 0 : void System::assemble_residual_derivatives (const ParameterVector &)
2554 : {
2555 0 : libmesh_not_implemented();
2556 : }
2557 :
2558 : inline
2559 0 : void System::disable_cache () { assemble_before_solve = true; }
2560 :
2561 : inline
2562 1702095 : unsigned int System::n_qois() const
2563 : {
2564 1702095 : libmesh_assert_equal_to(this->_qoi.size(), this->_qoi_error_estimates.size());
2565 :
2566 8194682 : return cast_int<unsigned int>(this->_qoi.size());
2567 : }
2568 :
2569 : inline
2570 : std::pair<unsigned int, Real>
2571 0 : System::sensitivity_solve (const ParameterVector &)
2572 : {
2573 0 : libmesh_not_implemented();
2574 : }
2575 :
2576 : inline
2577 : std::pair<unsigned int, Real>
2578 0 : System::weighted_sensitivity_solve (const ParameterVector &,
2579 : const ParameterVector &)
2580 : {
2581 0 : libmesh_not_implemented();
2582 : }
2583 :
2584 : inline
2585 : std::pair<unsigned int, Real>
2586 0 : System::adjoint_solve (const QoISet &)
2587 : {
2588 0 : libmesh_not_implemented();
2589 : }
2590 :
2591 : inline
2592 : std::pair<unsigned int, Real>
2593 0 : System::weighted_sensitivity_adjoint_solve (const ParameterVector &,
2594 : const ParameterVector &,
2595 : const QoISet &)
2596 : {
2597 0 : libmesh_not_implemented();
2598 : }
2599 :
2600 : inline
2601 : void
2602 0 : System::adjoint_qoi_parameter_sensitivity (const QoISet &,
2603 : const ParameterVector &,
2604 : SensitivityData &)
2605 : {
2606 0 : libmesh_not_implemented();
2607 : }
2608 :
2609 : inline
2610 : void
2611 0 : System::forward_qoi_parameter_sensitivity (const QoISet &,
2612 : const ParameterVector &,
2613 : SensitivityData &)
2614 : {
2615 0 : libmesh_not_implemented();
2616 : }
2617 :
2618 : inline
2619 : void
2620 0 : System::qoi_parameter_hessian(const QoISet &,
2621 : const ParameterVector &,
2622 : SensitivityData &)
2623 : {
2624 0 : libmesh_not_implemented();
2625 : }
2626 :
2627 : inline
2628 : void
2629 0 : System::qoi_parameter_hessian_vector_product(const QoISet &,
2630 : const ParameterVector &,
2631 : const ParameterVector &,
2632 : SensitivityData &)
2633 : {
2634 0 : libmesh_not_implemented();
2635 : }
2636 :
2637 : inline
2638 1052 : unsigned int System::n_matrices () const
2639 : {
2640 1052 : return cast_int<unsigned int>(_matrices.size());
2641 : }
2642 :
2643 : template <template <typename> class MatrixType>
2644 : inline
2645 : SparseMatrix<Number> &
2646 : System::add_matrix (std::string_view mat_name,
2647 : const ParallelType type)
2648 : {
2649 : // Return the matrix if it is already there.
2650 : auto it = this->_matrices.find(mat_name);
2651 : if (it != this->_matrices.end())
2652 : return *it->second;
2653 :
2654 : // Otherwise build the matrix to return.
2655 : auto pr = _matrices.emplace(mat_name, std::make_unique<MatrixType<Number>>(this->comm()));
2656 : _matrix_types.emplace(mat_name, type);
2657 :
2658 : SparseMatrix<Number> & mat = *(pr.first->second);
2659 :
2660 : // Initialize it first if we've already initialized the others.
2661 : this->late_matrix_init(mat, type);
2662 :
2663 : return mat;
2664 : }
2665 :
2666 : inline void
2667 : System::prefer_hash_table_matrix_assembly(const bool preference)
2668 : {
2669 : libmesh_error_msg_if(
2670 : _matrices_initialized,
2671 : "System::prefer_hash_table_matrix_assembly() should be called before matrices are initialized");
2672 : _prefer_hash_table_matrix_assembly = preference;
2673 : }
2674 : } // namespace libMesh
2675 :
2676 : #endif // LIBMESH_SYSTEM_H
|
