equation_systems.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2025 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 <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;

  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;

#ifdef LIBMESH_ENABLE_DEPRECATED

  void get_solution (std::vector<Number> & soln,
                     std::vector<std::string> & names) const;
#endif // LIBMESH_ENABLE_DEPRECATED

  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::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:
  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