static_condensation.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_STATIC_CONDENSATION_H
#define LIBMESH_STATIC_CONDENSATION_H

#include "libmesh/libmesh_config.h"

// subvectors currently only work with petsc and we rely on Eigen for local LU factorizations
#if defined(LIBMESH_HAVE_EIGEN) && defined(LIBMESH_HAVE_PETSC)

#include "libmesh/petsc_matrix_shell_matrix.h"
#include "libmesh/id_types.h"
#include "libmesh/libmesh_common.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/variable.h"
#include "libmesh/sparsity_pattern.h"

#include <unordered_map>
#include <memory>
#include <vector>

// Warnings about stack protection
#include "libmesh/ignore_warnings.h"
#include <Eigen/Dense>
#include "libmesh/restore_warnings.h"

namespace libMesh
{
class MeshBase;
class System;
class DofMap;
class Elem;
template <typename>
class LinearSolver;
template <typename>
class NumericVector;
template <typename>
class Preconditioner;
class StaticCondensationPreconditioner;
class StaticCondensationDofMap;

typedef Eigen::Matrix<Number, Eigen::Dynamic, Eigen::Dynamic> EigenMatrix;
typedef Eigen::Matrix<Number, Eigen::Dynamic, 1> EigenVector;

class StaticCondensation : public PetscMatrixShellMatrix<Number>
{
public:
  StaticCondensation(const MeshBase & mesh,
                     System & system,
                     const DofMap & full_dof_map,
                     StaticCondensationDofMap & reduced_dof_map);
  virtual ~StaticCondensation();

  //
  // SparseMatrix overrides
  //

  virtual SparseMatrix<Number> & operator=(const SparseMatrix<Number> &) override;

  virtual SolverPackage solver_package() override;

  virtual void init(const numeric_index_type m,
                    const numeric_index_type n,
                    const numeric_index_type m_l,
                    const numeric_index_type n_l,
                    const numeric_index_type nnz = 30,
                    const numeric_index_type noz = 10,
                    const numeric_index_type blocksize = 1) override;

  virtual void init(const ParallelType type) override;

  virtual bool initialized() const override;

  virtual void clear() noexcept override;

  virtual void zero() override;

  virtual std::unique_ptr<SparseMatrix<Number>> zero_clone() const override;

  virtual std::unique_ptr<SparseMatrix<Number>> clone() const override;

  virtual void close() override;

  virtual numeric_index_type m() const override;

  virtual numeric_index_type n() const override { return this->m(); }

  virtual numeric_index_type row_start() const override;

  virtual numeric_index_type row_stop() const override;

  virtual numeric_index_type col_start() const override { return this->row_start(); }

  virtual numeric_index_type col_stop() const override { return this->row_stop(); }

  virtual void
  set(const numeric_index_type i, const numeric_index_type j, const Number value) override;

  virtual void
  add(const numeric_index_type i, const numeric_index_type j, const Number value) override;

  virtual void add_matrix(const DenseMatrix<Number> & dm,
                          const std::vector<numeric_index_type> & rows,
                          const std::vector<numeric_index_type> & cols) override;

  virtual void add_matrix(const DenseMatrix<Number> & dm,
                          const std::vector<numeric_index_type> & dof_indices) override;

  virtual void add(const Number a, const SparseMatrix<Number> & X) override;

  virtual Number operator()(const numeric_index_type i, const numeric_index_type j) const override;

  virtual Real l1_norm() const override;

  virtual Real linfty_norm() const override;

  virtual bool closed() const override;

  virtual void print_personal(std::ostream & os = libMesh::out) const override;

  virtual void get_diagonal(NumericVector<Number> & dest) const override;

  virtual void get_transpose(SparseMatrix<Number> & dest) const override;

  virtual void get_row(numeric_index_type i,
                       std::vector<numeric_index_type> & indices,
                       std::vector<Number> & values) const override;

  //
  // Unique methods
  //

  void init();

  void setup();

  void apply(const NumericVector<Number> & full_rhs, NumericVector<Number> & full_sol);

  const SparseMatrix<Number> & get_condensed_mat() const;

  void set_current_elem(const Elem & elem);

  StaticCondensationPreconditioner & get_preconditioner() { return *_scp; }

  LinearSolver<Number> & reduced_system_solver();

  virtual bool require_sparsity_pattern() const override { return false; }

  void uncondensed_dofs_only() { _uncondensed_dofs_only = true; }

  void dont_condense_vars(const std::unordered_set<unsigned int> & vars);

private:
  static void set_local_vectors(const NumericVector<Number> & global_vector,
                                const std::vector<dof_id_type> & elem_dof_indices,
                                std::vector<Number> & elem_dof_values_vec,
                                EigenVector & elem_dof_values);

  void forward_elimination(const NumericVector<Number> & full_rhs);

  void backwards_substitution(const NumericVector<Number> & full_rhs,
                              NumericVector<Number> & full_sol);

  struct MatrixData
  {
    EigenMatrix Acc;
    EigenMatrix Acu;
    EigenMatrix Auc;
    EigenMatrix Auu;

    // Acc LU decompositions
    typename std::remove_const<decltype(Acc.partialPivLu())>::type AccFactor;
  };

  std::unordered_map<dof_id_type, MatrixData> _elem_to_matrix_data;

  const MeshBase & _mesh;
  System & _system;
  const DofMap & _full_dof_map;
  StaticCondensationDofMap & _reduced_dof_map;

  std::unique_ptr<SparseMatrix<Number>> _reduced_sys_mat;
  std::unique_ptr<NumericVector<Number>> _reduced_sol;
  std::unique_ptr<NumericVector<Number>> _reduced_rhs;
  std::unique_ptr<LinearSolver<Number>> _reduced_solver;
  // this is, in general, *not equal* to the system solution, e.g. this may correspond to the
  // solution for the Newton *update*
  std::unique_ptr<NumericVector<Number>> _ghosted_full_sol;

  dof_id_type _current_elem_id;

  DenseMatrix<Number> _size_one_mat;

  std::unique_ptr<StaticCondensationPreconditioner> _scp;

  bool _sc_is_initialized;

  bool _have_cached_values;

  ParallelType _parallel_type;

  bool _uncondensed_dofs_only;
};

inline const SparseMatrix<Number> & StaticCondensation::get_condensed_mat() const
{
  libmesh_assert(_reduced_sys_mat);
  return *_reduced_sys_mat;
}

inline LinearSolver<Number> & StaticCondensation::reduced_system_solver()
{
  libmesh_assert_msg(_reduced_solver, "Reduced system solver not built yet");
  return *_reduced_solver;
}

}

#else

#include "libmesh/sparse_matrix.h"
#include <unordered_set>

namespace libMesh
{
class MeshBase;
class System;
class DofMap;
class StaticCondensationPreconditioner;
class StaticCondensationDofMap;

class StaticCondensation : public SparseMatrix<Number>
{
public:
  StaticCondensation(const MeshBase &,
                     const System &,
                     const DofMap & full_dof_map,
                     StaticCondensationDofMap & reduced_dof_map);

  const std::unordered_set<unsigned int> & uncondensed_vars() const { libmesh_not_implemented(); }
  StaticCondensationPreconditioner & get_preconditioner() { libmesh_not_implemented(); }
  virtual SparseMatrix<Number> & operator=(const SparseMatrix<Number> &) override
  {
    libmesh_not_implemented();
  }
  virtual SolverPackage solver_package() override { libmesh_not_implemented(); }
  virtual void init(const numeric_index_type,
                    const numeric_index_type,
                    const numeric_index_type,
                    const numeric_index_type,
                    const numeric_index_type = 30,
                    const numeric_index_type = 10,
                    const numeric_index_type = 1) override
  {
    libmesh_not_implemented();
  }
  virtual void init(ParallelType) override { libmesh_not_implemented(); }
  virtual void clear() override { libmesh_not_implemented(); }
  virtual void zero() override { libmesh_not_implemented(); }
  virtual std::unique_ptr<SparseMatrix<Number>> zero_clone() const override
  {
    libmesh_not_implemented();
  }
  virtual std::unique_ptr<SparseMatrix<Number>> clone() const override
  {
    libmesh_not_implemented();
  }
  virtual void close() override { libmesh_not_implemented(); }
  virtual numeric_index_type m() const override { libmesh_not_implemented(); }
  virtual numeric_index_type n() const override { libmesh_not_implemented(); }
  virtual numeric_index_type row_start() const override { libmesh_not_implemented(); }
  virtual numeric_index_type row_stop() const override { libmesh_not_implemented(); }
  virtual numeric_index_type col_start() const override { libmesh_not_implemented(); }
  virtual numeric_index_type col_stop() const override { libmesh_not_implemented(); }
  virtual void set(const numeric_index_type, const numeric_index_type, const Number) override
  {
    libmesh_not_implemented();
  }
  virtual void add(const numeric_index_type, const numeric_index_type, const Number) override
  {
    libmesh_not_implemented();
  }
  virtual void add_matrix(const DenseMatrix<Number> &,
                          const std::vector<numeric_index_type> &,
                          const std::vector<numeric_index_type> &) override
  {
    libmesh_not_implemented();
  }
  virtual void add_matrix(const DenseMatrix<Number> &,
                          const std::vector<numeric_index_type> &) override
  {
    libmesh_not_implemented();
  }
  virtual void add(const Number, const SparseMatrix<Number> &) override
  {
    libmesh_not_implemented();
  }
  virtual Number operator()(const numeric_index_type, const numeric_index_type) const override
  {
    libmesh_not_implemented();
  }
  virtual Real l1_norm() const override { libmesh_not_implemented(); }
  virtual Real linfty_norm() const override { libmesh_not_implemented(); }
  virtual bool closed() const override { libmesh_not_implemented(); }
  virtual void print_personal(std::ostream & = libMesh::out) const override
  {
    libmesh_not_implemented();
  }
  virtual void get_diagonal(NumericVector<Number> &) const override { libmesh_not_implemented(); }
  virtual void get_transpose(SparseMatrix<Number> &) const override { libmesh_not_implemented(); }
  virtual void get_row(numeric_index_type,
                       std::vector<numeric_index_type> &,
                       std::vector<Number> &) const override
  {
    libmesh_not_implemented();
  }
  void init() { libmesh_not_implemented(); }
  void setup() { libmesh_not_implemented(); }
  void apply(const NumericVector<Number> &, NumericVector<Number> &) { libmesh_not_implemented(); }
};
}

#endif // LIBMESH_HAVE_EIGEN && LIBMESH_HAVE_PETSC
#endif // LIBMESH_STATIC_CONDENSATION_H