Discussion: Remove the inheritance of Epetra_Operator from SparseOperator?

src/adapter/4C_adapter_fld_fluid_fsi.cpp

@@ -511,7 +511,7 @@ void Adapter::FluidFSI::proj_vel_to_div_zero()
 
   std::shared_ptr<Core::LinAlg::Vector<double>> vmod =
       std::make_shared<Core::LinAlg::Vector<double>>(velnp()->get_map(), true);
-  B->Apply(*x, *vmod);
+  B->Apply(x->as_multi_vector(), vmod->as_multi_vector());
   write_access_velnp()->update(-1.0, *vmod, 1.0);
 }
 

src/ale/4C_ale.cpp

@@ -887,7 +887,7 @@ void ALE::AleLinear::evaluate_elements()
     validsysmat_ = true;
   }
   else if (system_matrix())
-    system_matrix()->Apply(*dispnp(), *write_access_residual());
+    system_matrix()->Apply(dispnp()->as_multi_vector(), write_access_residual()->as_multi_vector());
   else if (block_system_matrix())
     block_system_matrix()->Apply(*dispnp(), *write_access_residual());
   else

src/constraint/4C_constraint_manager.cpp

@@ -288,11 +288,7 @@ void Constraints::ConstrManager::evaluate_force_stiff(const double time,
   constrainterr_->update(scConMat, *referencevalues_, -1.0 * scConMat, *actvalues_, 0.0);
   actdisc_->clear_state();
   // finalize the constraint matrix
-  std::string label(constr_matrix_->Label());
-  if (label == "Core::LinAlg::BlockSparseMatrixBase")
-    constr_matrix_->complete();
-  else
-    constr_matrix_->complete(*constrmap_, *dofrowmap);
+  constr_matrix_->complete(*constrmap_, *dofrowmap);
 }
 
 /*----------------------------------------------------------------------*

src/contact/src/4C_contact_abstract_strategy.cpp

@@ -1482,7 +1482,8 @@ void CONTACT::AbstractStrategy::evaluate_relative_movement()
   for (int i = 0; i < (int)interfaces().size(); ++i) interfaces()[i]->assemble_slave_coord(xsmod);
 
   // in case of 3D dual quadratic case, slave coordinates xs are modified
-  if (is_dual_quad_slave_trafo()) invtrafo_->Apply(*xsmod, *xsmod);
+  if (is_dual_quad_slave_trafo())
+    invtrafo_->Apply(xsmod->as_multi_vector(), xsmod->as_multi_vector());
 
   // ATTENTION: for evaluate_relative_movement() we need the vector xsmod in
   // fully overlapping layout. Thus, export here. First, allreduce

src/core/linalg/src/sparse/4C_linalg_blocksparsematrix.hpp

@@ -27,7 +27,7 @@ namespace Core::LinAlg
     implemented in terms of the matrix blocks.
 
    */
-  class BlockSparseMatrixBase : public SparseOperator
+  class BlockSparseMatrixBase : public SparseOperator, public Epetra_Operator
   {
    public:
     /// constructor
@@ -121,6 +121,12 @@ namespace Core::LinAlg
     /// number of column blocks
     int cols() const { return domainmaps_.num_maps(); }
 
+    const Core::LinAlg::Map& range_map() const
+    {
+      FOUR_C_THROW("you need to specify which map you want. do not call this function.");
+      return *rangemaps_.map(0);
+    }
+
     /// range map for given row block
     const Core::LinAlg::Map& range_map(int r) const { return *rangemaps_.map(r); }
 
@@ -276,6 +282,9 @@ namespace Core::LinAlg
   class BlockSparseMatrix : public BlockSparseMatrixBase, public Strategy
   {
    public:
+    virtual Epetra_Operator& epetra_operator() { return this->epetra_operator(); };
+
+
     BlockSparseMatrix(const MultiMapExtractor& domainmaps, const MultiMapExtractor& rangemaps,
         int npr = 81, bool explicitdirichlet = true, bool savegraph = false);
 

src/core/linalg/src/sparse/4C_linalg_sparsematrix.cpp

@@ -1398,6 +1398,17 @@ int Core::LinAlg::SparseMatrix::Apply(const Epetra_MultiVector& X, Epetra_MultiV
   return sysmat_->Apply(X, Y);
 }
 
+int Core::LinAlg::SparseMatrix::Apply(
+    const Core::LinAlg::MultiVector<double>& X, Core::LinAlg::MultiVector<double>& Y) const
+{
+  return sysmat_->Apply(X.get_epetra_multi_vector(), Y.get_epetra_multi_vector());
+}
+
+int Core::LinAlg::SparseMatrix::Apply(
+    const Core::LinAlg::Vector<double>& X, Core::LinAlg::Vector<double>& Y) const
+{
+  return sysmat_->Apply(X.get_ref_of_epetra_vector(), Y.get_ref_of_epetra_vector());
+}
 
 /*----------------------------------------------------------------------*
  *----------------------------------------------------------------------*/

src/core/linalg/src/sparse/4C_linalg_sparsematrix.hpp

@@ -340,7 +340,7 @@ namespace Core::LinAlg
     //@{
 
     /// Returns a character string describing the operator.
-    const char* Label() const override;
+    const char* Label() const;
 
     //@}
 
@@ -375,7 +375,7 @@ namespace Core::LinAlg
 
     \note This method is here for performance reasons.
    */
-    Epetra_Operator& epetra_operator() override { return *sysmat_; }
+    Epetra_Operator& epetra_operator() { return *sysmat_; }
 
     /// return the internal Epetra matrix as Epetra_Operator
     const Epetra_Operator& epetra_operator() const { return *sysmat_; }
@@ -392,7 +392,7 @@ namespace Core::LinAlg
     //@{
 
     /// If set true, transpose of this operator will be applied.
-    int SetUseTranspose(bool UseTranspose) override;
+    int SetUseTranspose(bool UseTranspose);
 
     //@}
 
@@ -424,13 +424,18 @@ namespace Core::LinAlg
     //@{
 
     /// Returns the result of a Epetra_Operator applied to a Epetra_MultiVector X in Y.
-    int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const override;
+    int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
+
+    int Apply(
+        const Core::LinAlg::MultiVector<double>& X, Core::LinAlg::MultiVector<double>& Y) const;
+
+    int Apply(const Core::LinAlg::Vector<double>& X, Core::LinAlg::Vector<double>& Y) const;
 
     /// Returns the result of a Epetra_Operator inverse applied to an Epetra_MultiVector X in Y.
-    int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const override;
+    int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
 
     /// Returns the infinity norm of the global matrix.
-    double NormInf() const override;
+    double NormInf() const;
 
     /// Returns the one norm of the global matrix.
     double norm_one() const;
@@ -500,19 +505,19 @@ namespace Core::LinAlg
 
 
     /// Returns the current UseTranspose setting.
-    bool UseTranspose() const override;
+    bool UseTranspose() const;
 
     /// Returns true if the this object can provide an approximate Inf-norm, false otherwise.
-    bool HasNormInf() const override;
+    bool HasNormInf() const;
 
     /// Returns a pointer to the Epetra_Comm communicator associated with this operator.
-    const Epetra_Comm& Comm() const override;
+    const Epetra_Comm& Comm() const;
 
     /// Returns the Epetra_Map object associated with the domain of this operator.
-    const Epetra_Map& OperatorDomainMap() const override;
+    const Epetra_Map& OperatorDomainMap() const;
 
     /// Returns the Epetra_Map object associated with the range of this operator.
-    const Epetra_Map& OperatorRangeMap() const override;
+    const Epetra_Map& OperatorRangeMap() const;
 
     //@}
 

src/core/linalg/src/sparse/4C_linalg_sparseoperator.cpp

@@ -0,0 +1,17 @@
+// This file is part of 4C multiphysics licensed under the
+// GNU Lesser General Public License v3.0 or later.
+//
+// See the LICENSE.md file in the top-level for license information.
+//
+// SPDX-License-Identifier: LGPL-3.0-or-later
+
+#include "4C_linalg_sparseoperator.hpp"
+
+FOUR_C_NAMESPACE_OPEN
+
+namespace Core::LinAlg
+{
+  SparseOperator::~SparseOperator() {}
+}  // namespace Core::LinAlg
+
+FOUR_C_NAMESPACE_CLOSE

src/core/linalg/src/sparse/4C_linalg_sparseoperator.hpp

@@ -70,7 +70,7 @@ namespace Core::LinAlg
     is BlockSparseMatrix, a block matrix build from a list of SparseMatrix.
 
    */
-  class SparseOperator : public Epetra_Operator
+  class SparseOperator
   {
    public:
     /// return the internal Epetra_Operator
@@ -81,14 +81,17 @@ namespace Core::LinAlg
       \warning Only low level solver routines are interested in the internal
       Epetra_Operator.
      */
-    virtual Epetra_Operator& epetra_operator() { return *this; }
+    virtual ~SparseOperator();  // declaration only
+    virtual Epetra_Operator& epetra_operator() = 0;
 
     /// set matrix to zero
     virtual void zero() = 0;
 
     /// throw away the matrix and its graph and start anew
     virtual void reset() = 0;
 
+    virtual int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const = 0;
+
     /// Assemble a Core::LinAlg::SerialDenseMatrix into a matrix with striding
     /*!
 
@@ -211,6 +214,9 @@ namespace Core::LinAlg
     /// Returns the Epetra_Map object associated with the (full) domain of this operator.
     virtual const Map& domain_map() const = 0;
 
+    /// Returns the Epetra_Map object associated with the (full) range of this operator.
+    virtual const Map& range_map() const = 0;
+
     /// Add one operator to another
     virtual void add(const Core::LinAlg::SparseOperator& A, const bool transposeA,
         const double scalarA, const double scalarB) = 0;
@@ -237,4 +243,4 @@ namespace Core::LinAlg
 FOUR_C_NAMESPACE_CLOSE
 
 #endif
-/*LINALG_SPARSEOPERATOR_H_*/
+/*LINALG_SPARSEOPERATOR_H_*/

src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.cpp

@@ -48,8 +48,7 @@ void Core::LinearSolver::TekoPreconditioner::setup(Core::LinAlg::SparseOperator&
   auto xmlFileName = tekolist_.sublist("Teko Parameters").get<std::string>("TEKO_XML_FILE");
 
   Teuchos::ParameterList tekoParams;
-  auto comm = Core::Communication::to_teuchos_comm<int>(
-      Core::Communication::unpack_epetra_comm(matrix.Comm()));
+  auto comm = Core::Communication::to_teuchos_comm<int>(matrix.domain_map().get_comm());
   Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, Teuchos::Ptr(&tekoParams), *comm);
 
   auto A = std::dynamic_pointer_cast<Core::LinAlg::BlockSparseMatrixBase>(

src/fluid/4C_fluid_implicit_integration.cpp

@@ -1906,7 +1906,7 @@ void FLD::FluidImplicitTimeInt::evaluate_fluid_edge_based(
   if (systemmatrix1 != nullptr)
   {
     sysmat_linalg = std::make_shared<Core::LinAlg::SparseMatrix>(
-        systemmatrix1->OperatorRangeMap(), 256, true, false, Core::LinAlg::SparseMatrix::FE_MATRIX);
+        systemmatrix1->range_map(), 256, true, false, Core::LinAlg::SparseMatrix::FE_MATRIX);
   }
   else
     FOUR_C_THROW("sysmat is nullptr!");
@@ -2251,7 +2251,7 @@ void FLD::FluidImplicitTimeInt::check_matrix_nullspace()
 
     Core::LinAlg::Vector<double> result(c->get_map(), false);
 
-    sysmat_->Apply(c->get_epetra_multi_vector(), result);
+    sysmat_->Apply(c->get_epetra_multi_vector(), result.get_ref_of_epetra_vector());
 
     double norm = 1e9;
 
@@ -6344,8 +6344,8 @@ void FLD::FluidImplicitTimeInt::write_output_kinetic_energy()
 
   // compute kinetic energy
   double energy = 0.0;
-  Core::LinAlg::Vector<double> mtimesu(massmat_->OperatorRangeMap(), true);
-  massmat_->Apply(*velnp_, mtimesu);
+  Core::LinAlg::Vector<double> mtimesu(massmat_->range_map(), true);
+  massmat_->Apply(velnp_->get_ref_of_epetra_vector(), mtimesu.get_ref_of_epetra_vector());
   velnp_->dot(mtimesu, &energy);
   energy *= 0.5;
 

src/fsi/src/monolithic/nonlinear_solver/4C_fsi_nox_linearsystem.cpp

@@ -33,15 +33,16 @@ NOX::FSI::LinearSystem::LinearSystem(Teuchos::ParameterList& printParams,
     std::shared_ptr<Core::LinAlg::Solver> solver, const std::shared_ptr<NOX::Nln::Scaling> s)
     : utils_(printParams),
       jac_interface_ptr_(iJac),
-      jac_ptr_(J),
+      // jac_ptr_(Teuchos::rcp:_(J->epetra_operator())),
       operator_(J),
       scaling_(s),
       callcount_(0),
       solver_(solver),
       timer_("", true)
 {
   tmp_vector_ptr_ = std::make_shared<NOX::Nln::Vector>(cloneVector);
-
+  // TODO: Potential error here?
+  jac_ptr_ = std::shared_ptr<Epetra_Operator>(J, &J->epetra_operator());
   reset(linearSolverParams);
 }
 

src/solver_nonlin_nox/4C_solver_nonlin_nox_aux.cpp

@@ -64,12 +64,14 @@ NOX::Nln::LinSystem::OperatorType NOX::Nln::Aux::get_operator_type(
       const Core::LinAlg::BlockSparseMatrix<Core::LinAlg::DefaultBlockMatrixStrategy>*>(&op);
   if (testOperator != nullptr) return NOX::Nln::LinSystem::LinalgBlockSparseMatrix;
 
-  // Is it a LINALG_SparseMatrix?
-  testOperator = dynamic_cast<const Core::LinAlg::SparseMatrix*>(&op);
-  if (testOperator != nullptr) return NOX::Nln::LinSystem::LinalgSparseMatrix;
+  // TODO: Potential error here? - could probably fixed with usage of enum types.
+  //  Is it a LINALG_SparseMatrix
+  // testOperator = dynamic_cast<const Core::LinAlg::SparseMatrix*>(&op.epetra_operator());
+  // if (testOperator != nullptr)
+  return NOX::Nln::LinSystem::LinalgSparseMatrix;
 
   // Otherwise it must be a LINALG_SparseOperator
-  return NOX::Nln::LinSystem::LinalgSparseOperator;
+  // return NOX::Nln::LinSystem::LinalgSparseOperator;
 }
 
 /*----------------------------------------------------------------------------*

src/solver_nonlin_nox/4C_solver_nonlin_nox_linearsystem.cpp

@@ -195,7 +195,7 @@ bool NOX::Nln::LinearSystem::apply_jacobian_block(const NOX::Nln::Vector& input,
 bool NOX::Nln::LinearSystem::apply_jacobian(
     const NOX::Nln::Vector& input, NOX::Nln::Vector& result) const
 {
-  jacobian().SetUseTranspose(false);
+  jacobian().epetra_operator().SetUseTranspose(false);
   int status = jacobian().Apply(input.get_linalg_vector(), result.get_linalg_vector());
   return (status == 0);
 }
@@ -206,9 +206,9 @@ bool NOX::Nln::LinearSystem::apply_jacobian_transpose(
     const NOX::Nln::Vector& input, NOX::Nln::Vector& result) const
 {
   // Apply the Jacobian
-  jacobian().SetUseTranspose(true);
+  jacobian().epetra_operator().SetUseTranspose(true);
   int status = jacobian().Apply(input.get_linalg_vector(), result.get_linalg_vector());
-  jacobian().SetUseTranspose(false);
+  jacobian().epetra_operator().SetUseTranspose(false);
 
   return (status == 0);
 }
@@ -335,7 +335,8 @@ bool NOX::Nln::LinearSystem::apply_jacobian_inverse(Teuchos::ParameterList& line
 bool NOX::Nln::LinearSystem::compute_jacobian(const NOX::Nln::Vector& x)
 {
   prePostOperatorPtr_->run_pre_compute_jacobian(jacobian(), x.get_linalg_vector(), *this);
-  const bool success = jacInterfacePtr_->computeJacobian(x.get_linalg_vector(), jacobian());
+  const bool success =
+      jacInterfacePtr_->computeJacobian(x.get_linalg_vector(), jacobian().epetra_operator());
   prePostOperatorPtr_->run_post_compute_jacobian(jacobian(), x.get_linalg_vector(), *this);
 
   return success;
@@ -351,7 +352,8 @@ bool NOX::Nln::LinearSystem::compute_f_and_jacobian(
 
   const bool success =
       Teuchos::rcp_dynamic_cast<NOX::Nln::Interface::Jacobian>(jacInterfacePtr_, true)
-          ->compute_f_and_jacobian(x.get_linalg_vector(), rhs.get_linalg_vector(), jacobian());
+          ->compute_f_and_jacobian(
+              x.get_linalg_vector(), rhs.get_linalg_vector(), jacobian().epetra_operator());
 
   prePostOperatorPtr_->run_post_compute_f_and_jacobian(
       rhs.get_linalg_vector(), jacobian(), x.get_linalg_vector(), *this);
@@ -369,8 +371,8 @@ bool NOX::Nln::LinearSystem::compute_correction_system(const CorrectionType type
 
   const bool success =
       Teuchos::rcp_dynamic_cast<NOX::Nln::Interface::Jacobian>(jacInterfacePtr_, true)
-          ->compute_correction_system(
-              type, grp, x.get_linalg_vector(), rhs.get_linalg_vector(), jacobian());
+          ->compute_correction_system(type, grp, x.get_linalg_vector(), rhs.get_linalg_vector(),
+              jacobian().epetra_operator());
 
   prePostOperatorPtr_->run_post_compute_f_and_jacobian(
       rhs.get_linalg_vector(), jacobian(), x.get_linalg_vector(), *this);
@@ -478,14 +480,14 @@ NOX::Nln::LinearSystem::get_jacobian_interface() const
  *----------------------------------------------------------------------*/
 Teuchos::RCP<const Epetra_Operator> NOX::Nln::LinearSystem::get_jacobian_operator() const
 {
-  return jacobian_ptr();
+  return Teuchos::rcpFromRef(jac_ptr_->epetra_operator());
 }
 
 /*----------------------------------------------------------------------*
  *----------------------------------------------------------------------*/
 Teuchos::RCP<Epetra_Operator> NOX::Nln::LinearSystem::get_jacobian_operator()
 {
-  return jacobian_ptr();
+  return Teuchos::rcpFromRef(jac_ptr_->epetra_operator());
 }
 
 /*----------------------------------------------------------------------*
@@ -582,7 +584,7 @@ double NOX::Nln::LinearSystem::compute_serial_condition_number_of_jacobian(
     const LinSystem::ConditionNumber condnum_type) const
 {
   if (Core::Communication::num_mpi_ranks(
-          Core::Communication::unpack_epetra_comm(jacobian().Comm())) > 1)
+          Core::Communication::unpack_epetra_comm(jacobian().epetra_operator().Comm())) > 1)
     FOUR_C_THROW("Currently only one processor is supported!");
 
   Core::LinAlg::SerialDenseMatrix dense_jac;

src/solver_nonlin_nox/4C_solver_nonlin_nox_solver_ptc.cpp

@@ -9,6 +9,7 @@
 
 #include "4C_fem_geometry_intersection_math.hpp"
 #include "4C_linalg_sparsematrix.hpp"
+#include "4C_linalg_sparseoperator.hpp"
 #include "4C_linalg_utils_sparse_algebra_create.hpp"
 #include "4C_linalg_vector.hpp"
 #include "4C_solver_nonlin_nox_aux.hpp"

src/sti/4C_sti_monolithic.cpp

@@ -533,7 +533,7 @@ void STI::Monolithic::output_matrix_to_file(
   if (!sparseoperator->filled()) FOUR_C_THROW("Sparse operator must be filled for output!");
 
   // extract communicator
-  MPI_Comm comm = Core::Communication::unpack_epetra_comm(sparseoperator->Comm());
+  MPI_Comm comm = sparseoperator->range_map().get_comm();
 
   // determine whether sparse matrix or block sparse matrix should be output
   const auto sparsematrix =