system.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



#ifndef LIBMESH_SYSTEM_H
#define LIBMESH_SYSTEM_H

// Local Includes
#include "libmesh/elem_range.h"
#include "libmesh/enum_subset_solve_mode.h" // SUBSET_ZERO
#include "libmesh/enum_parallel_type.h" // PARALLEL
#include "libmesh/fem_function_base.h"
#include "libmesh/libmesh_common.h"
#include "libmesh/parallel_object.h"
#include "libmesh/parameters.h"
#include "libmesh/qoi_set.h"
#include "libmesh/reference_counted_object.h"
#include "libmesh/tensor_value.h" // For point_hessian
#include "libmesh/enum_matrix_build_type.h" // AUTOMATIC

// C++ includes
#include <cstddef>
#include <set>
#include <string>
#include <string_view>
#include <vector>
#include <memory>
#include <optional>

// This define may be useful in --disable-optional builds when it is
// possible that libmesh will not have any solvers available.
#if defined(LIBMESH_HAVE_PETSC) || \
  defined(LIBMESH_HAVE_EIGEN)   || \
  defined(LIBMESH_HAVE_LASPACK) || \
  defined(LIBMESH_TRILINOS_HAVE_AZTECOO)
#define LIBMESH_HAVE_SOLVER 1
#endif

namespace libMesh
{

// Forward Declarations
class System;
class EquationSystems;
class MeshBase;
class Xdr;
class DofMap;
template <typename Output> class FunctionBase;
class ParameterVector;
class Point;
class SensitivityData;
class Matrix;
template <typename T> class NumericVector;
template <typename T> class SparseMatrix;
template <typename T> class VectorValue;
typedef VectorValue<Number> NumberVectorValue;
typedef NumberVectorValue Gradient;
class SystemSubset;
class FEType;
class SystemNorm;
enum FEMNormType : int;
class Variable;
class VariableGroup;

class System : public ReferenceCountedObject<System>,
               public ParallelObject
{
public:

  System (EquationSystems & es,
          const std::string & name,
          const unsigned int number);

  System (const System &) = delete;
  System & operator= (const System &) = delete;
  System (System &&) = default;
  System & operator= (System &&) = delete;
  virtual ~System ();

  class Initialization
  {
  public:
    virtual ~Initialization () = default;

    virtual void initialize () = 0;
  };



  class Assembly
  {
  public:
    virtual ~Assembly () = default;

    virtual void assemble () = 0;
  };



  class Constraint
  {
  public:
    virtual ~Constraint () = default;

    virtual void constrain () = 0;
  };



  class QOI
  {
  public:
    virtual ~QOI () = default;

    virtual void qoi (const QoISet & qoi_indices) = 0;
  };



  class QOIDerivative
  {
  public:
    virtual ~QOIDerivative () = default;

    virtual void qoi_derivative (const QoISet & qoi_indices,
                                 bool include_liftfunc,
                                 bool apply_constraints) = 0;
  };

  typedef System sys_type;

  sys_type & system () { return *this; }

  virtual void clear ();

  void init ();

  virtual void reinit ();

  virtual void reinit_constraints ();

  virtual void reinit_mesh();

  bool is_initialized() const;

  virtual void update ();

  virtual void assemble ();

  virtual void assemble_qoi (const QoISet & qoi_indices = QoISet());

  virtual void assemble_qoi_derivative (const QoISet & qoi_indices = QoISet(),
                                        bool include_liftfunc = true,
                                        bool apply_constraints = true);

  virtual void assemble_residual_derivatives (const ParameterVector & parameters);

  virtual void restrict_solve_to (const SystemSubset * subset,
                                  const SubsetSolveMode subset_solve_mode=SUBSET_ZERO);

  virtual void solve () {}

  virtual std::pair<unsigned int, Real>
  sensitivity_solve (const ParameterVector & parameters);

  virtual std::pair<unsigned int, Real>
  weighted_sensitivity_solve (const ParameterVector & parameters,
                              const ParameterVector & weights);

  virtual std::pair<unsigned int, Real>
  adjoint_solve (const QoISet & qoi_indices = QoISet());

  virtual std::pair<unsigned int, Real>
  weighted_sensitivity_adjoint_solve (const ParameterVector & parameters,
                                      const ParameterVector & weights,
                                      const QoISet & qoi_indices = QoISet());
  bool is_adjoint_already_solved() const
  { return adjoint_already_solved;}

  void set_adjoint_already_solved(bool setting)
  { adjoint_already_solved = setting;}


  virtual void qoi_parameter_sensitivity (const QoISet & qoi_indices,
                                          const ParameterVector & parameters,
                                          SensitivityData & sensitivities);

  virtual void adjoint_qoi_parameter_sensitivity (const QoISet & qoi_indices,
                                                  const ParameterVector & parameters,
                                                  SensitivityData & sensitivities);

  virtual void forward_qoi_parameter_sensitivity (const QoISet & qoi_indices,
                                                  const ParameterVector & parameters,
                                                  SensitivityData & sensitivities);

  virtual void qoi_parameter_hessian(const QoISet & qoi_indices,
                                     const ParameterVector & parameters,
                                     SensitivityData & hessian);

  virtual void qoi_parameter_hessian_vector_product(const QoISet & qoi_indices,
                                                    const ParameterVector & parameters,
                                                    const ParameterVector & vector,
                                                    SensitivityData & product);

  virtual bool compare (const System & other_system,
                        const Real threshold,
                        const bool verbose) const;

  const std::string & name () const;

  virtual std::string system_type () const { return "Basic"; }

  void project_solution (FunctionBase<Number> * f,
                         FunctionBase<Gradient> * g = nullptr,
                         std::optional<ConstElemRange> active_local_range = std::nullopt,
                         std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;

  void project_solution (FEMFunctionBase<Number> * f,
                         FEMFunctionBase<Gradient> * g = nullptr,
                         std::optional<ConstElemRange> active_local_range = std::nullopt,
                         std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;

  typedef Number (*ValueFunctionPointer)(const Point & p,
                                         const Parameters & Parameters,
                                         const std::string & sys_name,
                                         const std::string & unknown_name);
  typedef Gradient (*GradientFunctionPointer)(const Point & p,
                                              const Parameters & parameters,
                                              const std::string & sys_name,
                                              const std::string & unknown_name);
  void project_solution (ValueFunctionPointer fptr,
                         GradientFunctionPointer gptr,
                         const Parameters & parameters,
                         std::optional<ConstElemRange> active_local_range = std::nullopt,
                         std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;

  void project_vector (NumericVector<Number> & new_vector,
                       FunctionBase<Number> * f,
                       FunctionBase<Gradient> * g = nullptr,
                       int is_adjoint = -1,
                       std::optional<ConstElemRange> active_local_range = std::nullopt,
                       std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;

  void project_vector (NumericVector<Number> & new_vector,
                       FEMFunctionBase<Number> * f,
                       FEMFunctionBase<Gradient> * g = nullptr,
                       int is_adjoint = -1,
                       std::optional<ConstElemRange> active_local_range = std::nullopt,
                       std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;

  void project_vector (ValueFunctionPointer fptr,
                       GradientFunctionPointer gptr,
                       const Parameters & parameters,
                       NumericVector<Number> & new_vector,
                       int is_adjoint = -1,
                       std::optional<ConstElemRange> active_local_range = std::nullopt,
                       std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;

  void boundary_project_solution (const std::set<boundary_id_type> & b,
                                  const std::vector<unsigned int> & variables,
                                  FunctionBase<Number> * f,
                                  FunctionBase<Gradient> * g = nullptr,
                                  std::optional<ConstElemRange> active_local_range = std::nullopt);

  void boundary_project_solution (const std::set<boundary_id_type> & b,
                                  const std::vector<unsigned int> & variables,
                                  ValueFunctionPointer fptr,
                                  GradientFunctionPointer gptr,
                                  const Parameters & parameters,
                                  std::optional<ConstElemRange> active_local_range = std::nullopt);

  void boundary_project_vector (const std::set<boundary_id_type> & b,
                                const std::vector<unsigned int> & variables,
                                NumericVector<Number> & new_vector,
                                FunctionBase<Number> * f,
                                FunctionBase<Gradient> * g = nullptr,
                                int is_adjoint = -1,
                                std::optional<ConstElemRange> active_local_range = std::nullopt) const;

  void boundary_project_vector (const std::set<boundary_id_type> & b,
                                const std::vector<unsigned int> & variables,
                                ValueFunctionPointer fptr,
                                GradientFunctionPointer gptr,
                                const Parameters & parameters,
                                NumericVector<Number> & new_vector,
                                int is_adjoint = -1,
                                std::optional<ConstElemRange> active_local_range = std::nullopt) const;

  unsigned int number () const;

  void update_global_solution (std::vector<Number> & global_soln) const;

  void update_global_solution (std::vector<Number> & global_soln,
                               const processor_id_type dest_proc) const;

  const MeshBase & get_mesh() const;

  MeshBase & get_mesh();

  const DofMap & get_dof_map() const;

  DofMap & get_dof_map();

  const EquationSystems & get_equation_systems() const { return _equation_systems; }

  EquationSystems & get_equation_systems() { return _equation_systems; }

  bool active () const;

  void activate ();

  void deactivate ();

  void set_basic_system_only ();

  typedef std::map<std::string, std::unique_ptr<NumericVector<Number>>, std::less<>>::iterator       vectors_iterator;
  typedef std::map<std::string, std::unique_ptr<NumericVector<Number>>, std::less<>>::const_iterator const_vectors_iterator;

  vectors_iterator vectors_begin ();

  const_vectors_iterator vectors_begin () const;

  vectors_iterator vectors_end ();

  const_vectors_iterator vectors_end () const;

  typedef std::map<std::string, std::unique_ptr<SparseMatrix<Number>>, std::less<>>::iterator        matrices_iterator;
  typedef std::map<std::string, std::unique_ptr<SparseMatrix<Number>>, std::less<>>::const_iterator  const_matrices_iterator;

  matrices_iterator matrices_begin ();

  const_matrices_iterator matrices_begin () const;

  matrices_iterator matrices_end ();

  const_matrices_iterator matrices_end () const;

  NumericVector<Number> & add_vector (std::string_view vec_name,
                                      const bool projections=true,
                                      const ParallelType type = PARALLEL);

  void remove_vector(std::string_view vec_name);

  bool & project_solution_on_reinit (void)
  { return _solution_projection; }

  bool have_vector (std::string_view vec_name) const;

  const NumericVector<Number> * request_vector (std::string_view vec_name) const;

  NumericVector<Number> * request_vector (std::string_view vec_name);

  const NumericVector<Number> * request_vector (const unsigned int vec_num) const;

  NumericVector<Number> * request_vector (const unsigned int vec_num);

  const NumericVector<Number> & get_vector (std::string_view vec_name) const;

  NumericVector<Number> & get_vector (std::string_view vec_name);

  const NumericVector<Number> & get_vector (const unsigned int vec_num) const;

  NumericVector<Number> & get_vector (const unsigned int vec_num);

  const std::string & vector_name (const unsigned int vec_num) const;

  const std::string & vector_name (const NumericVector<Number> & vec_reference) const;

  void set_vector_as_adjoint (const std::string & vec_name, int qoi_num);

  int vector_is_adjoint (std::string_view vec_name) const;

  void set_vector_preservation (const std::string & vec_name, bool preserve);

  bool vector_preservation (std::string_view vec_name) const;

  NumericVector<Number> & add_adjoint_solution(unsigned int i=0);

  NumericVector<Number> & get_adjoint_solution(unsigned int i=0);

  const NumericVector<Number> & get_adjoint_solution(unsigned int i=0) const;

  NumericVector<Number> & add_sensitivity_solution(unsigned int i=0);

  NumericVector<Number> & get_sensitivity_solution(unsigned int i=0);

  const NumericVector<Number> & get_sensitivity_solution(unsigned int i=0) const;

  NumericVector<Number> & add_weighted_sensitivity_adjoint_solution(unsigned int i=0);

  NumericVector<Number> & get_weighted_sensitivity_adjoint_solution(unsigned int i=0);

  const NumericVector<Number> & get_weighted_sensitivity_adjoint_solution(unsigned int i=0) const;

  NumericVector<Number> & add_weighted_sensitivity_solution();

  NumericVector<Number> & get_weighted_sensitivity_solution();

  const NumericVector<Number> & get_weighted_sensitivity_solution() const;

  NumericVector<Number> & add_adjoint_rhs(unsigned int i=0);

  NumericVector<Number> & get_adjoint_rhs(unsigned int i=0);

  const NumericVector<Number> & get_adjoint_rhs(unsigned int i=0) const;

  NumericVector<Number> & add_sensitivity_rhs(unsigned int i=0);

  NumericVector<Number> & get_sensitivity_rhs(unsigned int i=0);

  const NumericVector<Number> & get_sensitivity_rhs(unsigned int i=0) const;

  unsigned int n_vectors () const;

  unsigned int n_matrices () const;

  unsigned int n_vars() const;

  unsigned int n_variable_groups() const;

  unsigned int n_components() const;

  dof_id_type n_dofs() const;

  dof_id_type n_active_dofs() const;

  dof_id_type n_constrained_dofs() const;

  dof_id_type n_local_constrained_dofs() const;

  dof_id_type n_local_dofs() const;

  unsigned int add_variable (std::string_view var,
                             const FEType & type,
                             const std::set<subdomain_id_type> * const active_subdomains = nullptr);

  unsigned int add_variable (std::string_view var,
                             const Order order = FIRST,
                             const FEFamily = LAGRANGE,
                             const std::set<subdomain_id_type> * const active_subdomains = nullptr,
                             const bool p_refinement = true);

  unsigned int add_variables (const std::vector<std::string> & vars,
                              const FEType & type,
                              const std::set<subdomain_id_type> * const active_subdomains = nullptr);

  unsigned int add_variable_array (const std::vector<std::string> & vars,
                                   const FEType & type,
                                   const std::set<subdomain_id_type> * const active_subdomains = nullptr);

  unsigned int add_variables (const std::vector<std::string> & vars,
                              const Order order = FIRST,
                              const FEFamily = LAGRANGE,
                              const std::set<subdomain_id_type> * const active_subdomains = nullptr,
                              const bool p_refinement = true);

  const Variable & variable (unsigned int var) const;

  const VariableGroup & variable_group (unsigned int vg) const;

  bool has_variable(std::string_view var) const;

  const std::string & variable_name(const unsigned int i) const;

  unsigned int variable_number (std::string_view var) const;

  void get_all_variable_numbers(std::vector<unsigned int> & all_variable_numbers) const;

  unsigned int variable_scalar_number (std::string_view var,
                                       unsigned int component) const;

  unsigned int variable_scalar_number (unsigned int var_num,
                                       unsigned int component) const;


  const FEType & variable_type (const unsigned int i) const;

  const FEType & variable_type (std::string_view var) const;

  bool identify_variable_groups () const;

  void identify_variable_groups (const bool);

  Real calculate_norm(const NumericVector<Number> & v,
                      unsigned int var,
                      FEMNormType norm_type,
                      std::set<unsigned int> * skip_dimensions=nullptr) const;

  Real calculate_norm(const NumericVector<Number> & v,
                      const SystemNorm & norm,
                      std::set<unsigned int> * skip_dimensions=nullptr) const;

  void read_header (Xdr & io,
                    std::string_view version,
                    const bool read_header=true,
                    const bool read_additional_data=true,
                    const bool read_legacy_format=false);

  template <typename ValType>
  void read_serialized_data (Xdr & io,
                             const bool read_additional_data=true);
  void read_serialized_data (Xdr & io,
                             const bool read_additional_data=true)
  { read_serialized_data<Number>(io, read_additional_data); }

  template <typename InValType>
  std::size_t read_serialized_vectors (Xdr & io,
                                       const std::vector<NumericVector<Number> *> & vectors) const;

  std::size_t read_serialized_vectors (Xdr & io,
                                       const std::vector<NumericVector<Number> *> & vectors) const
  { return read_serialized_vectors<Number>(io, vectors); }

  template <typename InValType>
  void read_parallel_data (Xdr & io,
                           const bool read_additional_data);

  void read_parallel_data (Xdr & io,
                           const bool read_additional_data)
  { read_parallel_data<Number>(io, read_additional_data); }

  void write_header (Xdr & io,
                     std::string_view version,
                     const bool write_additional_data) const;

  void write_serialized_data (Xdr & io,
                              const bool write_additional_data = true) const;

  std::size_t write_serialized_vectors (Xdr & io,
                                        const std::vector<const NumericVector<Number> *> & vectors) const;

  void write_parallel_data (Xdr & io,
                            const bool write_additional_data) const;

  std::string get_info () const;

  void attach_init_function (void fptr(EquationSystems & es,
                                       const std::string & name));

  void attach_init_object (Initialization & init);

  void attach_assemble_function (void fptr(EquationSystems & es,
                                           const std::string & name));

  void attach_assemble_object (Assembly & assemble);

  void attach_constraint_function (void fptr(EquationSystems & es,
                                             const std::string & name));

  void attach_constraint_object (Constraint & constrain);

  bool has_constraint_object () const;

  Constraint& get_constraint_object ();

  void attach_QOI_function (void fptr(EquationSystems & es,
                                      const std::string & name,
                                      const QoISet & qoi_indices));

  void attach_QOI_object (QOI & qoi);

  void attach_QOI_derivative (void fptr(EquationSystems & es,
                                        const std::string & name,
                                        const QoISet & qoi_indices,
                                        bool include_liftfunc,
                                        bool apply_constraints));

  void attach_QOI_derivative_object (QOIDerivative & qoi_derivative);

  virtual void user_initialization ();

  virtual void user_assembly ();

  virtual void user_constrain ();

  virtual void user_QOI (const QoISet & qoi_indices);

  virtual void user_QOI_derivative (const QoISet & qoi_indices = QoISet(),
                                    bool include_liftfunc = true,
                                    bool apply_constraints = true);

  virtual void re_update ();

  virtual void restrict_vectors ();

  virtual void prolong_vectors ();

  Parameters parameters;

  bool assemble_before_solve;

  virtual void disable_cache ();

  bool use_fixed_solution;

  int extra_quadrature_order;


  //--------------------------------------------------
  // The solution and solution access members

  Number current_solution (const dof_id_type global_dof_number) const;

  std::unique_ptr<NumericVector<Number>> solution;

  std::unique_ptr<NumericVector<Number>> current_local_solution;

  Real time;

  unsigned int n_qois() const;

public:

  void init_qois(unsigned int n_qois);

  void set_qoi(unsigned int qoi_index, Number qoi_value);
  Number get_qoi_value(unsigned int qoi_index) const;

  void set_qoi(std::vector<Number> new_qoi);

  std::vector<Number> get_qoi_values() const;

  void set_qoi_error_estimate(unsigned int qoi_index, Number qoi_error_estimate);
  Number get_qoi_error_estimate_value(unsigned int qoi_index) const;

  Number point_value(unsigned int var, const Point & p,
                     const bool insist_on_success = true,
                     const NumericVector<Number> * sol = nullptr) const;

  Number point_value(unsigned int var, const Point & p, const Elem & e,
                     const NumericVector<Number> * sol = nullptr) const;

  Number point_value(unsigned int var, const Point & p, const Elem * e) const;

  Number point_value(unsigned int var, const Point & p, const NumericVector<Number> * sol) const;

  Gradient point_gradient(unsigned int var, const Point & p,
                          const bool insist_on_success = true,
                          const NumericVector<Number> * sol = nullptr) const;

  Gradient point_gradient(unsigned int var, const Point & p, const Elem & e,
                          const NumericVector<Number> * sol = nullptr) const;

  Gradient point_gradient(unsigned int var, const Point & p, const Elem * e) const;

  Gradient point_gradient(unsigned int var, const Point & p, const NumericVector<Number> * sol) const;

  Tensor point_hessian(unsigned int var, const Point & p,
                       const bool insist_on_success = true,
                       const NumericVector<Number> * sol = nullptr) const;

  Tensor point_hessian(unsigned int var, const Point & p, const Elem & e,
                       const NumericVector<Number> * sol = nullptr) const;

  Tensor point_hessian(unsigned int var, const Point & p, const Elem * e) const;

  Tensor point_hessian(unsigned int var, const Point & p, const NumericVector<Number> * sol) const;


  void local_dof_indices (const unsigned int var,
                          std::set<dof_id_type> & var_indices) const;

  void zero_variable (NumericVector<Number> & v, unsigned int var_num) const;

  bool get_project_with_constraints()
  {
    return project_with_constraints;
  }

  void set_project_with_constraints(bool _project_with_constraints)
  {
    project_with_constraints = _project_with_constraints;
  }

  bool & hide_output() { return _hide_output; }

#ifdef LIBMESH_HAVE_METAPHYSICL

  void projection_matrix (SparseMatrix<Number> & proj_mat) const;
#endif // LIBMESH_HAVE_METAPHYSICL

  SparseMatrix<Number> & add_matrix (std::string_view mat_name,
                                     ParallelType type = PARALLEL,
                                     MatrixBuildType mat_build_type = MatrixBuildType::AUTOMATIC);

  template <template <typename> class>
  SparseMatrix<Number> & add_matrix (std::string_view mat_name, ParallelType = PARALLEL);

  SparseMatrix<Number> & add_matrix (std::string_view mat_name,
                                     std::unique_ptr<SparseMatrix<Number>> matrix,
                                     ParallelType type = PARALLEL);

  void remove_matrix(std::string_view mat_name);

  inline bool have_matrix (std::string_view mat_name) const { return _matrices.count(mat_name); }

  const SparseMatrix<Number> * request_matrix (std::string_view mat_name) const;

  SparseMatrix<Number> * request_matrix (std::string_view mat_name);

  const SparseMatrix<Number> & get_matrix (std::string_view mat_name) const;

  SparseMatrix<Number> & get_matrix (std::string_view mat_name);

  void prefer_hash_table_matrix_assembly(bool preference);

  void prefix_with_name(bool value) { _prefix_with_name = value; }

  [[nodiscard]] bool prefix_with_name() const { return _prefix_with_name; }

  std::string prefix() const { return this->name() + "_"; }

   virtual void create_static_condensation();

   bool has_static_condensation() const;

  /*
   * If we have e.g. a element space constrained by spline values, we
   * can directly project functions only on the constrained basis; to
   * get consistent constraining values we have to solve for them.
   *
   * Constrain the new vector using the requested adjoint rather than
   * primal constraints if is_adjoint is non-negative.
   */
  void solve_for_unconstrained_dofs (NumericVector<Number> &,
                                     int is_adjoint = -1) const;

protected:

  virtual void init_data ();

  virtual void add_matrices() {}

  virtual void init_matrices ();

  bool can_add_matrices() const { return !_matrices_initialized; }

  void project_vector (NumericVector<Number> &,
                       int is_adjoint = -1,
                       std::optional<ConstElemRange> active_local_range = std::nullopt,
                       std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;

  void project_vector (const NumericVector<Number> &,
                       NumericVector<Number> &,
                       int is_adjoint = -1,
                       std::optional<ConstElemRange> active_local_range = std::nullopt,
                       std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;

  virtual bool condense_constrained_dofs() const { return false; }

private:
  void late_matrix_init(SparseMatrix<Number> & mat,
                        ParallelType type);

  Real discrete_var_norm (const NumericVector<Number> & v,
                          unsigned int var,
                          FEMNormType norm_type) const;

  template <typename iterator_type, typename InValType>
  std::size_t read_serialized_blocked_dof_objects (const dof_id_type n_objects,
                                                   const iterator_type begin,
                                                   const iterator_type end,
                                                   const InValType dummy,
                                                   Xdr & io,
                                                   const std::vector<NumericVector<Number> *> & vecs,
                                                   const unsigned int var_to_read=libMesh::invalid_uint) const;

  unsigned int read_SCALAR_dofs (const unsigned int var,
                                 Xdr & io,
                                 NumericVector<Number> * vec) const;

  template <typename InValType>
  numeric_index_type read_serialized_vector (Xdr & io,
                                             NumericVector<Number> * vec);

  numeric_index_type read_serialized_vector (Xdr & io,
                                             NumericVector<Number> & vec)
  { return read_serialized_vector<Number>(io, &vec); }

  template <typename iterator_type>
  std::size_t write_serialized_blocked_dof_objects (const std::vector<const NumericVector<Number> *> & vecs,
                                                    const dof_id_type n_objects,
                                                    const iterator_type begin,
                                                    const iterator_type end,
                                                    Xdr & io,
                                                    const unsigned int var_to_write=libMesh::invalid_uint) const;

  unsigned int write_SCALAR_dofs (const NumericVector<Number> & vec,
                                  const unsigned int var,
                                  Xdr & io) const;

  dof_id_type write_serialized_vector (Xdr & io,
                                       const NumericVector<Number> & vec) const;

  void (* _init_system_function) (EquationSystems & es,
                                  const std::string & name);

  Initialization * _init_system_object;

  void (* _assemble_system_function) (EquationSystems & es,
                                      const std::string & name);

  Assembly * _assemble_system_object;

  void (* _constrain_system_function) (EquationSystems & es,
                                       const std::string & name);

  Constraint * _constrain_system_object;

  void (* _qoi_evaluate_function) (EquationSystems & es,
                                   const std::string & name,
                                   const QoISet & qoi_indices);

  QOI * _qoi_evaluate_object;

  void (* _qoi_evaluate_derivative_function) (EquationSystems & es,
                                              const std::string & name,
                                              const QoISet & qoi_indices,
                                              bool include_liftfunc,
                                              bool apply_constraints);

  QOIDerivative * _qoi_evaluate_derivative_object;

  std::unique_ptr<DofMap> _dof_map;

  EquationSystems & _equation_systems;

  MeshBase & _mesh;

  const std::string _sys_name;

  const unsigned int _sys_number;

  bool _active;

  std::map<std::string, std::unique_ptr<NumericVector<Number>>, std::less<>> _vectors;

  std::map<std::string, bool, std::less<>> _vector_projections;

  std::map<std::string, int, std::less<>> _vector_is_adjoint;

  std::map<std::string, std::unique_ptr<SparseMatrix<Number>>, std::less<>> _matrices;

  std::map<std::string, ParallelType, std::less<>> _matrix_types;

  bool _matrices_initialized;

  bool _solution_projection;

  bool _basic_system_only;

  bool _is_initialized;

  unsigned int _additional_data_written;

  std::vector<unsigned int> _written_var_indices;

  bool adjoint_already_solved;

  bool _hide_output;

  bool project_with_constraints;

  bool _prefer_hash_table_matrix_assembly;

  bool _require_sparsity_pattern;

  bool _prefix_with_name;

  std::vector<Number> _qoi;

  std::vector<Number> _qoi_error_estimates;
};



// ------------------------------------------------------------
// System inline methods
inline
const std::string & System::name() const
{
  return _sys_name;
}



inline
unsigned int System::number() const
{
  return _sys_number;
}



inline
const MeshBase & System::get_mesh() const
{
  return _mesh;
}



inline
MeshBase & System::get_mesh()
{
  return _mesh;
}



inline
const DofMap & System::get_dof_map() const
{
  return *_dof_map;
}



inline
DofMap & System::get_dof_map()
{
  return *_dof_map;
}



inline
bool System::active() const
{
  return _active;
}



inline
void System::activate ()
{
  _active = true;
}



inline
void System::deactivate ()
{
  _active = false;
}



inline
bool System::is_initialized () const
{
  return _is_initialized;
}



inline
void System::set_basic_system_only ()
{
  _basic_system_only = true;
}



inline
unsigned int
System::variable_scalar_number (std::string_view var,
                                unsigned int component) const
{
  return variable_scalar_number(this->variable_number(var), component);
}



inline
dof_id_type System::n_active_dofs() const
{
  return this->n_dofs() - this->n_constrained_dofs();
}



inline
bool System::have_vector (std::string_view vec_name) const
{
  return (_vectors.count(vec_name));
}



inline
unsigned int System::n_vectors () const
{
  return cast_int<unsigned int>(_vectors.size());
}

inline
System::vectors_iterator System::vectors_begin ()
{
  return _vectors.begin();
}

inline
System::const_vectors_iterator System::vectors_begin () const
{
  return _vectors.begin();
}

inline
System::vectors_iterator System::vectors_end ()
{
  return _vectors.end();
}

inline
System::const_vectors_iterator System::vectors_end () const
{
  return _vectors.end();
}

inline
System::matrices_iterator System::matrices_begin ()
{
  return _matrices.begin();
}

inline
System::const_matrices_iterator System::matrices_begin () const
{
  return _matrices.begin();
}

inline
System::matrices_iterator System::matrices_end ()
{
  return _matrices.end();
}

inline
System::const_matrices_iterator System::matrices_end () const
{
  return _matrices.end();
}

inline
void System::assemble_residual_derivatives (const ParameterVector &)
{
  libmesh_not_implemented();
}

inline
void System::disable_cache () { assemble_before_solve = true; }

inline
unsigned int System::n_qois() const
{
  libmesh_assert_equal_to(this->_qoi.size(), this->_qoi_error_estimates.size());

  return cast_int<unsigned int>(this->_qoi.size());
}

inline
std::pair<unsigned int, Real>
System::sensitivity_solve (const ParameterVector &)
{
  libmesh_not_implemented();
}

inline
std::pair<unsigned int, Real>
System::weighted_sensitivity_solve (const ParameterVector &,
                                    const ParameterVector &)
{
  libmesh_not_implemented();
}

inline
std::pair<unsigned int, Real>
System::adjoint_solve (const QoISet &)
{
  libmesh_not_implemented();
}

inline
std::pair<unsigned int, Real>
System::weighted_sensitivity_adjoint_solve (const ParameterVector &,
                                            const ParameterVector &,
                                            const QoISet &)
{
  libmesh_not_implemented();
}

inline
void
System::adjoint_qoi_parameter_sensitivity (const QoISet &,
                                           const ParameterVector &,
                                           SensitivityData &)
{
  libmesh_not_implemented();
}

inline
void
System::forward_qoi_parameter_sensitivity (const QoISet &,
                                           const ParameterVector &,
                                           SensitivityData &)
{
  libmesh_not_implemented();
}

inline
void
System::qoi_parameter_hessian(const QoISet &,
                              const ParameterVector &,
                              SensitivityData &)
{
  libmesh_not_implemented();
}

inline
void
System::qoi_parameter_hessian_vector_product(const QoISet &,
                                             const ParameterVector &,
                                             const ParameterVector &,
                                             SensitivityData &)
{
  libmesh_not_implemented();
}

inline
unsigned int System::n_matrices () const
{
  return cast_int<unsigned int>(_matrices.size());
}

template <template <typename> class MatrixType>
inline
SparseMatrix<Number> &
System::add_matrix (std::string_view mat_name,
                    const ParallelType type)
{
  // Return the matrix if it is already there.
  auto it = this->_matrices.find(mat_name);
  if (it != this->_matrices.end())
    return *it->second;

  // Otherwise build the matrix to return.
  auto pr = _matrices.emplace(mat_name, std::make_unique<MatrixType<Number>>(this->comm()));
  _matrix_types.emplace(mat_name, type);

  SparseMatrix<Number> & mat = *(pr.first->second);

  // Initialize it first if we've already initialized the others.
  this->late_matrix_init(mat, type);

  return mat;
}

inline void
System::prefer_hash_table_matrix_assembly(const bool preference)
{
  libmesh_error_msg_if(
      _matrices_initialized,
      "System::prefer_hash_table_matrix_assembly() should be called before matrices are initialized");
  _prefer_hash_table_matrix_assembly = preference;
}
} // namespace libMesh

#endif // LIBMESH_SYSTEM_H