Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2025 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             : #ifndef LIBMESH_CONDENSED_EIGEN_SYSTEM_H
      20             : #define LIBMESH_CONDENSED_EIGEN_SYSTEM_H
      21             : 
      22             : #include "libmesh/libmesh_config.h"
      23             : 
      24             : // Currently, the EigenSystem should only be available
      25             : // if SLEPc support is enabled.
      26             : #if defined(LIBMESH_HAVE_SLEPC)
      27             : 
      28             : // Local Includes
      29             : #include "libmesh/eigen_system.h"
      30             : #include "libmesh/sparse_matrix.h"
      31             : 
      32             : namespace libMesh
      33             : {
      34             : 
      35             : /**
      36             :  * This class extends EigenSystem to allow a simple way of solving
      37             :  * (standard or generalized) eigenvalue problems in the case where
      38             :  * we want to remove certain degrees of freedom from the system.
      39             :  * This is useful, for example, in the case that one wants to solve
      40             :  * eigenvalue problems with Dirichlet boundary conditions.
      41             :  *
      42             :  * \author David Knezevic
      43             :  * \date 2011
      44             :  * \brief Extends EigenSystem to allow certain DOFs to be condensed out.
      45             :  */
      46          20 : class CondensedEigenSystem : public EigenSystem
      47             : {
      48             : public:
      49             : 
      50             :   /**
      51             :    * Constructor.
      52             :    */
      53             :   CondensedEigenSystem (EquationSystems & es,
      54             :                         const std::string & name_in,
      55             :                         const unsigned int number_in);
      56             : 
      57             :   /**
      58             :    * Special functions.
      59             :    * - This class has the same restrictions/defaults as its base class.
      60             :    * - The destructor is defaulted out-of-line.
      61             :    */
      62             :   CondensedEigenSystem (const CondensedEigenSystem &) = delete;
      63             :   CondensedEigenSystem & operator= (const CondensedEigenSystem &) = delete;
      64             :   CondensedEigenSystem (CondensedEigenSystem &&) = default;
      65             :   CondensedEigenSystem & operator= (CondensedEigenSystem &&) = delete;
      66             :   virtual ~CondensedEigenSystem ();
      67             : 
      68             :   /**
      69             :    * The type of system.
      70             :    */
      71             :   typedef CondensedEigenSystem sys_type;
      72             : 
      73             :   /**
      74             :    * The type of the parent
      75             :    */
      76             :   typedef EigenSystem Parent;
      77             : 
      78             :   /**
      79             :    * \returns A reference to *this.
      80             :    */
      81             :   sys_type & system () { return *this; }
      82             : 
      83             :   /**
      84             :    * Loop over the dofs on each processor to initialize the list
      85             :    * of non-condensed dofs. These are the dofs in the system that
      86             :    * are not contained in \p global_dirichlet_dofs_set and are not
      87             :    * subject to constraints due to adaptive mesh hanging nodes,
      88             :    * periodic boundary conditions, or Dirichlet boundary conditions.
      89             :    *
      90             :    * Most users will not need to use the \p global_condensed_dofs_set
      91             :    * argument; simply call initialize_condensed_dofs() after any time
      92             :    * the EquationSystems (and therefore its constraint equations) gets
      93             :    * initialized or reinitialized.
      94             :    */
      95             :   void initialize_condensed_dofs(const std::set<dof_id_type> &
      96             :                                  global_condensed_dofs_set =
      97             :                                  std::set<dof_id_type>());
      98             : 
      99             :   /**
     100             :    * \returns The global number of non-condensed dofs in the system.
     101             :    */
     102             :   dof_id_type n_global_non_condensed_dofs() const;
     103             : 
     104             :   /**
     105             :    * Override to solve the condensed eigenproblem with
     106             :    * the dofs in local_non_condensed_dofs_vector
     107             :    * stripped out of the system matrices on each processor.
     108             :    */
     109             :   virtual void solve() override;
     110             : 
     111             :   virtual void clear () override;
     112             : 
     113             :   /**
     114             :    * Override \p get_eigenpair() to retrieve the eigenpair for
     115             :    * the condensed eigensolve. We only set the non-condensed
     116             :    * entries of the solution vector (the condensed
     117             :    * entries are set to zero by default).
     118             :    */
     119             :   virtual std::pair<Real, Real> get_eigenpair(dof_id_type i) override;
     120             : 
     121           0 :   bool has_condensed_matrix_A() const { return _condensed_matrix_A.get(); }
     122           0 :   bool has_condensed_matrix_B() const { return _condensed_matrix_B.get(); }
     123           0 :   bool has_condensed_precond_matrix() const { return _condensed_precond_matrix.get(); }
     124             : 
     125             :   /**
     126             :    * \returns The system matrix used for standard eigenvalue problems
     127             :    */
     128             :   SparseMatrix<Number> & get_condensed_matrix_A();
     129             : 
     130             :   /**
     131             :    * \returns The second system matrix used for generalized eigenvalue problems.
     132             :    */
     133             :   SparseMatrix<Number> & get_condensed_matrix_B();
     134             : 
     135             :   /**
     136             :    * \returns The condensed preconditioning matrix
     137             :    */
     138             :   SparseMatrix<Number> & get_condensed_precond_matrix();
     139             : 
     140             :   /**
     141             :    * Copy a logically sub-vector into a super-vector. The sub-vector local size should correspond to
     142             :    * the size of our \p local_non_condensed_dofs_vector, e.g. the sub-vector should represent
     143             :    * non-condensed degree of freedom. The \p super should contain both condensed and non-condensed
     144             :    * dofs
     145             :    */
     146             :   void copy_sub_to_super(const NumericVector<Number> & sub, NumericVector<Number> & super);
     147             : 
     148             :   /**
     149             :    * Copy a logically super-vector into a sub-vector. The \p super vector should represent all dofs,
     150             :    * both condensed and non-condensed. The \p sub should contain only non-condensed dofs, e.g. its
     151             :    * local size should match the size of our \p local_non_condensed_dofs_vector
     152             :    */
     153             :   void copy_super_to_sub(NumericVector<Number> & super, NumericVector<Number> & sub);
     154             : 
     155             :   /**
     156             :    * Copy a logically super-matrix into a sub-matrix. The \p super matrix should represent all dofs,
     157             :    * both condensed and non-condensed. The \p sub should contain only non-condensed dofs, e.g. its
     158             :    * local row ownership size should match the size of our \p local_non_condensed_dofs_vector
     159             :    */
     160             :   void copy_super_to_sub(const SparseMatrix<Number> & super, SparseMatrix<Number> & sub);
     161             : 
     162             :   /**
     163             :    * Instructs not to create the condensed submatrices from the global matrices right before the
     164             :    * solve. This API should be used when assembly occurs during callbacks during the solve, e.g. for
     165             :    * a nonlinear eigen problem, solved with a Newton-based solver to determine the fundamental mode,
     166             :    * in which function and matrix evalations are liable to change with every nonlinear iteration.
     167             :    * Moreover, calling \p create_submatrix wipes away callbacks associated with the condensed matrix
     168             :    */
     169             :   void dont_create_submatrices_in_solve() { _create_submatrices_in_solve = false; }
     170             : 
     171             : #ifdef LIBMESH_ENABLE_DEPRECATED
     172             :   /**
     173             :    * The (condensed) system matrix for standard eigenvalue problems.
     174             :    *
     175             :    * Public access to this member variable will be deprecated in
     176             :    * the future! Use get_condensed_matrix_A() instead.
     177             :    */
     178             :   SparseMatrix<Number> * condensed_matrix_A;
     179             : 
     180             :   /**
     181             :    * A second (condensed) system matrix for generalized eigenvalue problems.
     182             :    *
     183             :    * Public access to this member variable will be deprecated in
     184             :    * the future! Use get_condensed_matrix_B() instead.
     185             :    */
     186             :   SparseMatrix<Number> * condensed_matrix_B;
     187             : #endif
     188             : 
     189             :   /**
     190             :    * Vector storing the local dof indices that will not be condensed.
     191             :    * All dofs that are not in this vector will be eliminated from
     192             :    * the system when we perform a solve.
     193             :    */
     194             :   std::vector<dof_id_type> local_non_condensed_dofs_vector;
     195             : 
     196             :   /**
     197             :    * Initializes the condensed matrices. This API should be used if the condensed matrices will be
     198             :    * assembled during the solve as opposed to pre-solve
     199             :    */
     200             :   void initialize_condensed_matrices();
     201             : 
     202             :   /**
     203             :    * @returns Whether there are condensed degrees of freedom
     204             :    */
     205             :   bool have_condensed_dofs() const
     206             :   { libmesh_assert(_condensed_dofs_initialized); return _have_condensed_dofs; }
     207             : 
     208             :   virtual void reinit() override;
     209             : 
     210             : protected:
     211             :   virtual void add_matrices () override;
     212             : 
     213             : private:
     214         208 :   virtual bool condense_constrained_dofs() const override final { return true; }
     215             : 
     216             : #ifdef LIBMESH_ENABLE_DEPRECATED
     217             :   void set_raw_pointers();
     218             : #endif
     219             : 
     220             :   /**
     221             :    * A private flag to indicate whether the condensed dofs
     222             :    * have been initialized.
     223             :    */
     224             :   bool _condensed_dofs_initialized;
     225             : 
     226             :   /**
     227             :    * Whether there are any condensed degrees of freedom
     228             :    */
     229             :   bool _have_condensed_dofs;
     230             : 
     231             :   /**
     232             :    * Denotes whether to create the condensed submatrices from the global matrices in the solve
     233             :    */
     234             :   bool _create_submatrices_in_solve;
     235             : 
     236             :   /**
     237             :    * The (condensed) system matrix for standard eigenvalue problems.
     238             :    */
     239             :   std::unique_ptr<SparseMatrix<Number>> _condensed_matrix_A;
     240             : 
     241             :   /**
     242             :    * A second (condensed) system matrix for generalized eigenvalue problems.
     243             :    */
     244             :   std::unique_ptr<SparseMatrix<Number>> _condensed_matrix_B;
     245             : 
     246             :   /**
     247             :    * The condensed preconditioning matrix
     248             :    */
     249             :   std::unique_ptr<SparseMatrix<Number>> _condensed_precond_matrix;
     250             : };
     251             : 
     252             : 
     253             : } // namespace libMesh
     254             : 
     255             : 
     256             : #endif // LIBMESH_HAVE_SLEPC
     257             : 
     258             : #endif // LIBMESH_CONDENSED_EIGEN_SYSTEM_H