/* SPDX-FileCopyrightText: Copyright (c) 2022 - 2026, the adamantine authors.
 * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 */

#ifndef MECHANICAL_OPERATOR_HH
#define MECHANICAL_OPERATOR_HH

#include <BodyForce.hh>
#include <MaterialProperty.hh>

#include <deal.II/base/memory_space.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/hp/q_collection.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/la_parallel_vector.h>
#if DEAL_II_VERSION_GTE(9, 7, 0) && defined(DEAL_II_TRILINOS_WITH_TPETRA)
#include <deal.II/lac/trilinos_tpetra_precondition.h>
#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
#else
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#endif

#include <boost/property_tree/ptree.hpp>

namespace adamantine
{
/**
 * This class is the operator associated with the solid mechanics equations.
 * The class is templated on the MemorySpace because it use MaterialProperty
 * which itself is templated on the MemorySpace but the operator is CPU only.
 */
template <int dim, int n_materials, int p_order, typename MaterialStates,
          typename MemorySpaceType>
class MechanicalOperator
{
public:
#if DEAL_II_VERSION_GTE(9, 7, 0) && defined(DEAL_II_TRILINOS_WITH_TPETRA)
  using TrilinosMatrixType =
      dealii::LinearAlgebra::TpetraWrappers::SparseMatrix<
          double, dealii::MemorySpace::Default>;
  using TrilinosPreconditionerType =
      dealii::LinearAlgebra::TpetraWrappers::PreconditionAMGMueLu<
          double, dealii::MemorySpace::Default>;
#else
  using TrilinosMatrixType = dealii::TrilinosWrappers::SparseMatrix;
  using TrilinosPreconditionerType = dealii::TrilinosWrappers::PreconditionAMG;
#endif

  /**
   * Constructor. If the initial temperature is negative, the simulation is
   * mechanical only. Otherwise, we solve a thermo-mechanical problem.
   */
  MechanicalOperator(MPI_Comm const &communicator,
                     MaterialProperty<dim, n_materials, p_order, MaterialStates,
                                      MemorySpaceType> &material_properties,
                     std::vector<double> const &reference_temperatures);

  void reinit(dealii::DoFHandler<dim> const &dof_handler,
              dealii::AffineConstraints<double> const &affine_constraints,
              dealii::hp::QCollection<dim> const &quad,
              std::vector<std::shared_ptr<BodyForce<dim>>> const &body_forces =
                  std::vector<std::shared_ptr<BodyForce<dim>>>());

  /**
   * Update the DoFHandler used by ThermalPhysics and update the temperature.
   */
  void update_temperature(
      dealii::DoFHandler<dim> const &thermal_dof_handler,
      dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host> const
          &temperature,
      std::vector<bool> const &has_melted);

  /**
   * Assemble the right-hand-side.
   */
  void
  assemble_rhs(std::vector<std::shared_ptr<BodyForce<dim>>> const &body_forces);

  /**
   * Return the right-hand-side.
   */
  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host> const &
  rhs() const;

  /**
   * Return the matrix of the system.
   */
  TrilinosMatrixType const &system_matrix() const;

  /**
   * Return the preconditioner for the system matrix.
   */
  TrilinosPreconditionerType const &preconditioner() const;

private:
  /**
   * Assemble the sparse matrix of the system.
   * @note The 2D case does not represent any physical model but it is
   * convenient for testing.
   */
  void assemble_matrix();

  /**
   * MPI communicator.
   */
  MPI_Comm const &_communicator;
  /**
   * List of initial temperatures of the material. If the length of the vector
   * is nonzero, we solve a thermo-mechanical problem.
   */
  std::vector<double> _reference_temperatures;
  /**
   * Reference to the MaterialProperty from MechanicalPhysics.
   */
  MaterialProperty<dim, n_materials, p_order, MaterialStates, MemorySpaceType>
      &_material_properties;
  /**
   * Non-owning pointer to the DoFHandler from MechanicalPhysics
   */
  dealii::DoFHandler<dim> const *_dof_handler = nullptr;
  /**
   * Non-owning pointer to the DoFHandler from ThermalPhysics
   */
  dealii::DoFHandler<dim> const *_thermal_dof_handler = nullptr;
  /**
   * Non-owning pointer to the AffineConstraints from MechanicalPhysics
   */
  dealii::AffineConstraints<double> const *_affine_constraints = nullptr;
  /**
   * Non-owning pointer to the QCollection from MechanicalPhysics
   */
  dealii::hp::QCollection<dim> const *_q_collection = nullptr;
  /**
   * Right-hand-side of the mechanical problem.
   */
  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>
      _system_rhs;
  /**
   * Matrix of the mechanical problem.
   */
  TrilinosMatrixType _system_matrix;
  /**
   * Preconditioner associated with the system matrix
   */
  TrilinosPreconditionerType _preconditioner;
  /**
   * Temperature of the material.
   */
  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>
      _temperature;
  /**
   * Indicator variable for whether a point has ever been above the solidus. The
   * value is false for material that has not yet melted and true for material
   * that has melted.
   */
  std::vector<bool> _has_melted;
};

template <int dim, int n_materials, int p_order, typename MaterialStates,
          typename MemorySpaceType>
inline dealii::LA::distributed::Vector<double,
                                       dealii::MemorySpace::Host> const &
MechanicalOperator<dim, n_materials, p_order, MaterialStates,
                   MemorySpaceType>::rhs() const
{
  return _system_rhs;
}

template <int dim, int n_materials, int p_order, typename MaterialStates,
          typename MemorySpaceType>
inline typename MechanicalOperator<dim, n_materials, p_order, MaterialStates,
                                   MemorySpaceType>::TrilinosMatrixType const &
MechanicalOperator<dim, n_materials, p_order, MaterialStates,
                   MemorySpaceType>::system_matrix() const
{
  return _system_matrix;
}

template <int dim, int n_materials, int p_order, typename MaterialStates,
          typename MemorySpaceType>
inline typename MechanicalOperator<
    dim, n_materials, p_order, MaterialStates,
    MemorySpaceType>::TrilinosPreconditionerType const &
MechanicalOperator<dim, n_materials, p_order, MaterialStates,
                   MemorySpaceType>::preconditioner() const
{
  return _preconditioner;
}
} // namespace adamantine
#endif