equation_systems.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_EQUATION_SYSTEMS_H
#define LIBMESH_EQUATION_SYSTEMS_H
 
// Local Includes
#include "libmesh/libmesh_common.h"
#include "libmesh/parameters.h"
#include "libmesh/system.h"
#include "libmesh/parallel_object.h"
 
// HP aCC needs these for some reason
#ifdef __HP_aCC
# include "libmesh/frequency_system.h"
# include "libmesh/transient_system.h"
# include "libmesh/newmark_system.h"
# include "libmesh/steady_system.h"
#endif
 
// C++ includes
#include <cstddef>
#include <functional>
#include <map>
#include <set>
#include <string>
#include <string_view>
#include <vector>
#include <memory>
 
namespace libMesh
{
 
// Forward Declarations
class Elem;
class MeshBase;
enum XdrMODE : int;
 
class EquationSystems : public ReferenceCountedObject<EquationSystems>,
                        public ParallelObject
 
{
public:
 
  enum ReadFlags { READ_HEADER           = 1,
                   READ_DATA             = 2,
                   READ_ADDITIONAL_DATA  = 4,
                   READ_LEGACY_FORMAT    = 8,
                   TRY_READ_IFEMS        = 16,
                   READ_BASIC_ONLY       = 32 };
 
  enum WriteFlags { WRITE_DATA             = 1,
                    WRITE_ADDITIONAL_DATA  = 2,
                    WRITE_PARALLEL_FILES   = 4,
                    WRITE_SERIAL_FILES     = 8 };
 
  EquationSystems (MeshBase & mesh);
 
  virtual ~EquationSystems ();
 
  virtual void clear ();
 
  virtual void init ();
 
  virtual void reinit ();
 
  virtual void reinit_mesh ();
 
  virtual void enable_default_ghosting (bool enable);
 
  void update ();
 
  unsigned int n_systems() const;
 
  bool has_system (std::string_view name) const;
 
  template <typename T_sys>
  const T_sys & get_system (std::string_view name) const;
 
  template <typename T_sys>
  T_sys & get_system (std::string_view name);
 
  template <typename T_sys>
  const T_sys & get_system (const unsigned int num) const;
 
  template <typename T_sys>
  T_sys & get_system (const unsigned int num);
 
  const System & get_system (std::string_view name) const;
 
  System & get_system (std::string_view name);
 
  const System & get_system (const unsigned int num) const;
 
  System & get_system (const unsigned int num);
 
  virtual System & add_system (std::string_view system_type,
                               std::string_view name);
 
  template <typename T_sys>
  T_sys & add_system (std::string_view name);
 
  unsigned int n_vars () const;
 
  std::size_t n_dofs () const;
 
  std::size_t n_active_dofs() const;
 
  virtual void solve ();
 
  virtual void adjoint_solve (const QoISet & qoi_indices = QoISet());
 
  virtual void sensitivity_solve (const ParameterVector & parameters);
 
  void build_variable_names (std::vector<std::string> & var_names,
                             const FEType * type=nullptr,
                             const std::set<std::string> * system_names=nullptr) const;
 
  static bool is_elemental_data_fe_type (const FEType & type);
 
  void build_elemental_data_variable_names
    (std::vector<std::string> & var_names,
     const std::set<std::string> * system_names=nullptr) const;
 
  void build_solution_vector (std::vector<Number> & soln,
                              std::string_view system_name,
                              std::string_view variable_name = "all_vars") const;
 
  void build_solution_vector (std::vector<Number> & soln,
                              const std::set<std::string> * system_names=nullptr,
                              bool add_sides=false) const;
 
  std::unique_ptr<NumericVector<Number>>
  build_parallel_solution_vector(const std::set<std::string> * system_names=nullptr,
                                 bool add_sides=false) const;
 
  void get_vars_active_subdomains(const std::vector<std::string> & names,
                                  std::vector<std::set<subdomain_id_type>> & vars_active_subdomains) const;
 
  void build_elemental_solution_vector (std::vector<Number> & soln,
                                        std::vector<std::string> & names) const;
 
  std::vector<std::pair<unsigned int, unsigned int>>
  find_variable_numbers (std::vector<std::string> & names,
                         const FEType * type=nullptr,
                         const std::vector<FEType> * types=nullptr) const;
 
  std::vector<std::pair<unsigned int, unsigned int>>
  find_elemental_data_variable_numbers (std::vector<std::string> & names) const;
 
  std::unique_ptr<NumericVector<Number>>
  build_parallel_elemental_solution_vector (std::vector<std::string> & names) const;
 
  void build_discontinuous_solution_vector
  (std::vector<Number> & soln,
   const std::set<std::string> * system_names = nullptr,
   const std::vector<std::string> * var_names = nullptr,
   bool vertices_only = false,
   bool add_sides = false) const;
 
  /*
   * Returns true iff the given side of the given element is *never*
   * added to output from that element, because it is considered to be
   * redundant with respect to the same data added from the
   * neighboring element sharing that side.  This helper function is
   * used with the add_sides option when building solution vectors
   * here and when outputting solution vectors in MeshOutput
   * (currently just Exodus) I/O.
   */
  static bool redundant_added_side(const Elem & elem, unsigned int side);
 
 
  template <typename InValType = Number>
  void read (std::string_view name,
             const XdrMODE,
             const unsigned int read_flags=(READ_HEADER | READ_DATA),
             bool partition_agnostic = true);
 
  template <typename InValType = Number>
  void read (std::string_view name,
             const unsigned int read_flags=(READ_HEADER | READ_DATA),
             bool partition_agnostic = true);
 
  template <typename InValType = Number>
  void read (Xdr & io,
             std::function<std::unique_ptr<Xdr>()> & local_io_functor,
             const unsigned int read_flags=(READ_HEADER | READ_DATA),
             bool partition_agnostic = true);
 
  void write (std::string_view name,
              const XdrMODE,
              const unsigned int write_flags=(WRITE_DATA),
              bool partition_agnostic = true) const;
 
  void write (std::string_view name,
              const unsigned int write_flags=(WRITE_DATA),
              bool partition_agnostic = true) const;
 
  void write (std::ostream name,
              const unsigned int write_flags=(WRITE_DATA),
              bool partition_agnostic = true) const;
 
  void write (Xdr & io,
              const unsigned int write_flags=(WRITE_DATA),
              bool partition_agnostic = true,
              Xdr * const local_io = nullptr) const;
 
  virtual bool compare (const EquationSystems & other_es,
                        const Real threshold,
                        const bool verbose) const;
 
  virtual std::string get_info() const;
 
  void print_info (std::ostream & os=libMesh::out) const;
 
  friend std::ostream & operator << (std::ostream & os,
                                     const EquationSystems & es);
 
  const MeshBase & get_mesh() const;
 
  MeshBase & get_mesh();
 
  void allgather ();
 
  void enable_refine_in_reinit() { this->_refine_in_reinit = true; }
 
  void disable_refine_in_reinit() { this->_refine_in_reinit = false; }
 
  bool refine_in_reinit_flag() { return this->_refine_in_reinit; }
 
  bool reinit_solutions ();
 
  virtual void reinit_systems ();
 
  Parameters parameters;
 
 
protected:
 
 
  MeshBase & _mesh;
 
  std::map<std::string, std::unique_ptr<System>, std::less<>> _systems;
 
  bool _refine_in_reinit;
 
  bool _enable_default_ghosting;
 
private:
  std::vector<std::pair<unsigned int, unsigned int>>
  find_variable_numbers_by_predicate
    (std::vector<std::string> & names,
     const std::function<bool(const FEType &)> & type_filter) const;
 
  void _add_system_to_nodes_and_elems();
 
  void _remove_default_ghosting(unsigned int sys_num);
};
 
 
 
// ------------------------------------------------------------
// EquationSystems inline methods
inline
const MeshBase & EquationSystems::get_mesh () const
{
  return _mesh;
}
 
 
 
inline
MeshBase & EquationSystems::get_mesh ()
{
  return _mesh;
}
 
 
inline
unsigned int EquationSystems::n_systems () const
{
  return cast_int<unsigned int>(_systems.size());
}
 
 
 
 
template <typename T_sys>
inline
T_sys & EquationSystems::add_system (std::string_view name)
{
  if (!_systems.count(name))
    {
      const unsigned int sys_num = this->n_systems();
 
      auto result = _systems.emplace
        (name, std::make_unique<T_sys>(*this, std::string(name),
                                       sys_num));
 
      if (!_enable_default_ghosting)
        this->_remove_default_ghosting(sys_num);
 
      // Tell all the \p DofObject entities to add a system.
      this->_add_system_to_nodes_and_elems();
 
      // Return reference to newly added item
      auto it = result.first;
      auto & sys_ptr = it->second;
      return cast_ref<T_sys &>(*sys_ptr);
    }
  else
    {
      // We now allow redundant add_system calls, to make it
      // easier to load data from files for user-derived system
      // subclasses
      return this->get_system<T_sys>(name);
    }
}
 
 
 
inline
bool EquationSystems::has_system (std::string_view name) const
{
  if (_systems.find(name) == _systems.end())
    return false;
  return true;
}
 
 
 
 
template <typename T_sys>
inline
const T_sys & EquationSystems::get_system (const unsigned int num) const
{
  libmesh_assert_less (num, this->n_systems());
 
  for (auto & pr : _systems)
    {
      const auto & sys_ptr = pr.second;
      if (sys_ptr->number() == num)
        return cast_ref<const T_sys &>(*sys_ptr);
    }
  // Error if we made it here
  libmesh_error_msg("ERROR: no system number " << num << " found!");
}
 
 
 
 
template <typename T_sys>
inline
T_sys & EquationSystems::get_system (const unsigned int num)
{
  libmesh_assert_less (num, this->n_systems());
 
  for (auto & pr : _systems)
    {
      auto & sys_ptr = pr.second;
      if (sys_ptr->number() == num)
        return cast_ref<T_sys &>(*sys_ptr);
    }
 
  // Error if we made it here
  libmesh_error_msg("ERROR: no system number " << num << " found!");
}
 
 
 
 
 
 
template <typename T_sys>
inline
const T_sys & EquationSystems::get_system (std::string_view name) const
{
  auto pos = _systems.find(name);
 
  // Check for errors
  libmesh_error_msg_if(pos == _systems.end(), "ERROR: no system named \"" << name << "\" found!");
 
  // Attempt dynamic cast
  const auto & sys_ptr = pos->second;
  return cast_ref<const T_sys &>(*sys_ptr);
}
 
 
 
 
 
 
template <typename T_sys>
inline
T_sys & EquationSystems::get_system (std::string_view name)
{
  auto pos = _systems.find(name);
 
  // Check for errors
  libmesh_error_msg_if(pos == _systems.end(), "ERROR: no system named " << name << " found!");
 
  // Attempt dynamic cast
  auto & sys_ptr = pos->second;
  return cast_ref<T_sys &>(*sys_ptr);
}
 
 
 
 
 
 
 
inline
const System & EquationSystems::get_system (std::string_view name) const
{
  return this->get_system<System>(name);
}
 
 
 
inline
System & EquationSystems::get_system (std::string_view name)
{
  return this->get_system<System>(name);
}
 
 
 
inline
const System & EquationSystems::get_system (const unsigned int num) const
{
  return this->get_system<System>(num);
}
 
 
 
inline
System & EquationSystems::get_system (const unsigned int num)
{
  return this->get_system<System>(num);
}
 
 
} // namespace libMesh
 
 
#endif // LIBMESH_EQUATION_SYSTEMS_H