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             : #ifndef LIBMESH_REDUCED_SPARSE_MATRIX_H
      19             : #define LIBMESH_REDUCED_SPARSE_MATRIX_H
      20             : 
      21             : #include "libmesh/libmesh_config.h"
      22             : 
      23             : #ifdef LIBMESH_HAVE_PETSC
      24             : 
      25             : #include "libmesh/petsc_matrix_base.h"
      26             : 
      27             : #include <memory>
      28             : #include <vector>
      29             : 
      30             : namespace libMesh
      31             : {
      32             : class ReducedDofMap;
      33             : template <typename T>
      34             : class PetscMatrix;
      35             : 
      36             : /**
      37             :  * A matrix on the reduced space of a \p ReducedDofMap, assembled with full degree of freedom
      38             :  * indices.
      39             :  *
      40             :  * Insertion (\p set(), \p add(), \p add_matrix()) takes indices of the full \p DofMap, because
      41             :  * element assembly produces them. Each index is mapped to its reduced index, and entries in rows
      42             :  * or columns that the \p ReducedDofMap does not retain are dropped. Dropping is exact when the
      43             :  * element contributions already carry the reduction, e.g. \f$ C^T K C \f$ for eliminated
      44             :  * constraints. Every other operation uses reduced indices: sizes, row and column ranges, entry
      45             :  * access, and the PETSc \p Mat returned by \p mat(), which is what solvers see.
      46             :  */
      47             : template <typename T>
      48          60 : class ReducedSparseMatrix : public PetscMatrixBase<T>
      49             : {
      50             : public:
      51             :   ReducedSparseMatrix(const Parallel::Communicator & comm_in, const ReducedDofMap & reduced_dof_map);
      52             : 
      53             :   virtual ~ReducedSparseMatrix();
      54             : 
      55             :   ReducedSparseMatrix(ReducedSparseMatrix &&) = delete;
      56             :   ReducedSparseMatrix(const ReducedSparseMatrix &) = delete;
      57             :   ReducedSparseMatrix & operator=(ReducedSparseMatrix &&) = delete;
      58             :   ReducedSparseMatrix & operator=(const ReducedSparseMatrix &) = delete;
      59             : 
      60             :   virtual SparseMatrix<T> & operator=(const SparseMatrix<T> & v) override;
      61             : 
      62             :   /**
      63             :    * Allocates the reduced matrix from the \p ReducedDofMap sparsity pattern
      64             :    */
      65             :   virtual void init(ParallelType type = PARALLEL) override;
      66             : 
      67             :   /**
      68             :    * Arguments are reduced sizes
      69             :    */
      70             :   virtual void init(const numeric_index_type m,
      71             :                     const numeric_index_type n,
      72             :                     const numeric_index_type m_l,
      73             :                     const numeric_index_type n_l,
      74             :                     const numeric_index_type nnz = 30,
      75             :                     const numeric_index_type noz = 10,
      76             :                     const numeric_index_type blocksize = 1) override;
      77             : 
      78          36 :   virtual bool require_sparsity_pattern() const override { return false; }
      79             : 
      80             :   virtual void clear() noexcept override;
      81             : 
      82             :   virtual void zero() override;
      83             : 
      84             :   virtual std::unique_ptr<SparseMatrix<T>> zero_clone() const override;
      85             : 
      86             :   virtual std::unique_ptr<SparseMatrix<T>> clone() const override;
      87             : 
      88             :   virtual void close() override;
      89             : 
      90             :   virtual void flush() override;
      91             : 
      92             :   /**
      93             :    * \p i and \p j are full indices
      94             :    */
      95             :   virtual void set(const numeric_index_type i, const numeric_index_type j, const T value) override;
      96             : 
      97             :   /**
      98             :    * \p i and \p j are full indices
      99             :    */
     100             :   virtual void add(const numeric_index_type i, const numeric_index_type j, const T value) override;
     101             : 
     102             :   /**
     103             :    * \p rows and \p cols are full indices
     104             :    */
     105             :   virtual void add_matrix(const DenseMatrix<T> & dm,
     106             :                           const std::vector<numeric_index_type> & rows,
     107             :                           const std::vector<numeric_index_type> & cols) override;
     108             : 
     109             :   /**
     110             :    * \p dof_indices are full indices
     111             :    */
     112             :   virtual void add_matrix(const DenseMatrix<T> & dm,
     113             :                           const std::vector<numeric_index_type> & dof_indices) override;
     114             : 
     115             :   virtual void add(const T a, const SparseMatrix<T> & X) override;
     116             : 
     117             :   /**
     118             :    * \p i and \p j are reduced indices
     119             :    */
     120             :   virtual T operator()(const numeric_index_type i, const numeric_index_type j) const override;
     121             : 
     122             :   virtual Real l1_norm() const override;
     123             : 
     124             :   virtual Real linfty_norm() const override;
     125             : 
     126             :   virtual void print_personal(std::ostream & os = libMesh::out) const override;
     127             : 
     128             :   /**
     129             :    * Fills the reduced-size vector \p dest with the diagonal
     130             :    */
     131             :   virtual void get_diagonal(NumericVector<T> & dest) const override;
     132             : 
     133             :   virtual void get_transpose(SparseMatrix<T> & dest) const override;
     134             : 
     135             :   /**
     136             :    * \p i and the returned \p indices are reduced indices
     137             :    */
     138             :   virtual void get_row(numeric_index_type i,
     139             :                        std::vector<numeric_index_type> & indices,
     140             :                        std::vector<T> & values) const override;
     141             : 
     142             :   virtual void scale(const T scale) override;
     143             : 
     144             :   /**
     145             :    * \returns The reduced degree of freedom map
     146             :    */
     147           0 :   const ReducedDofMap & reduced_dof_map() const { return _reduced_dof_map; }
     148             : 
     149             :   /**
     150             :    * \returns The wrapped matrix on the reduced space
     151             :    */
     152             :   PetscMatrix<T> & reduced_matrix();
     153             : 
     154             : private:
     155             :   /**
     156             :    * Maps the full indices \p full to reduced indices in \p reduced, with \p -1 for unretained
     157             :    * ones, which PETSc ignores on insertion
     158             :    */
     159             :   void to_reduced(const std::vector<numeric_index_type> & full,
     160             :                   std::vector<PetscInt> & reduced) const;
     161             : 
     162             :   /**
     163             :    * Wraps the reduced matrix's \p Mat and sets it as this matrix's context
     164             :    */
     165             :   void adopt_reduced_matrix();
     166             : 
     167             :   const ReducedDofMap & _reduced_dof_map;
     168             : 
     169             :   /// The matrix on the reduced space
     170             :   std::unique_ptr<PetscMatrix<T>> _reduced;
     171             : 
     172             :   /// Work vectors of reduced row and column indices
     173             :   mutable std::vector<PetscInt> _rows, _cols;
     174             : };
     175             : 
     176             : } // namespace libMesh
     177             : 
     178             : #endif // LIBMESH_HAVE_PETSC
     179             : #endif // LIBMESH_REDUCED_SPARSE_MATRIX_H