// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.

#ifndef MFEML2_MORTAR_INTEGRATOR_HPP
#define MFEML2_MORTAR_INTEGRATOR_HPP

#include "../fem.hpp"

namespace mfem
{

/*!
 * @brief Interface for mortar element assembly.
 * The MortarIntegrator interface is used for performing Petrov-Galerkin
 * finite element assembly on intersections between elements.
 * The quadrature rules are to be generated by a cut algorithm (e.g.,
 * mfem::Cut). The quadrature rules are defined in the respective trial and test
 * reference frames. Trial and test spaces can be associated with different
 * element shapes (e.g., triangles and quadrilaterals) and different polynomial
 * orders (e.g., 1 and 4). This class is designed to work in conjunction with
 * the MFEM/moonolith module but it can be used also for other applications.
 */
class MortarIntegrator
{
public:
   /*!
    * @brief Implements the assembly routine
    * @param trial is the master/source element
    * @param trial_ir is the quadrature formula for evaluating quantities within
    * the trial element
    * @param trial_Trans the geometric transformation of the trial element
    * @param test  is the slave/destination element
    * @param test_ir is the quadrature formula for evaluating quantities within
    * the test element
    * @param test_Trans the geometric transformation of the test element
    * @param elemmat the result of the assembly
    */
   virtual void AssembleElementMatrix(const FiniteElement &trial,
                                      const IntegrationRule &trial_ir,
                                      ElementTransformation &trial_Trans,
                                      const FiniteElement &test,
                                      const IntegrationRule &test_ir,
                                      ElementTransformation &test_Trans,
                                      DenseMatrix &elemmat) = 0;

   /*!
    * @return the additional orders of quadrature required by the integrator.
    * It is 0 by default, override method to change that.
    */
   virtual int GetQuadratureOrder() const { return 0; }

   virtual ~MortarIntegrator() {}

   /*!
    * @return true if it uses vector fe and false otherwise
    */
   virtual bool is_vector_fe() const = 0;
};

/*!
 * @brief Integrator for scalar finite elements
 * \f$ (u, v)_{L^2(\mathcal{T}_m \cap \mathcal{T}_s)}, u \in U(\mathcal{T}_m )
 * and v \in V(\mathcal{T}_s ) \f$
 */
class L2MortarIntegrator : public MortarIntegrator
{
public:
   void AssembleElementMatrix(const FiniteElement &trial,
                              const IntegrationRule &trial_ir,
                              ElementTransformation &trial_Trans,
                              const FiniteElement &test,
                              const IntegrationRule &test_ir,
                              ElementTransformation &test_Trans,
                              DenseMatrix &elemmat) override;

   inline bool is_vector_fe() const override { return false; }
};

/*!
 * @brief Integrator for vector finite elements. Experimental.
 * \f$ (u, v)_{L^2(\mathcal{T}_m \cap \mathcal{T}_s)}, u \in U(\mathcal{T}_m )
 * and v \in V(\mathcal{T}_s ) \f$
 */
class VectorL2MortarIntegrator : public MortarIntegrator
{
public:
#ifndef MFEM_THREAD_SAFE
   Vector shape;
   Vector D;
   DenseMatrix K;
   DenseMatrix test_vshape;
   DenseMatrix trial_vshape;
#endif

public:
   VectorL2MortarIntegrator() { Init(NULL, NULL, NULL); }
   VectorL2MortarIntegrator(Coefficient *_q) { Init(_q, NULL, NULL); }
   VectorL2MortarIntegrator(VectorCoefficient *_vq) { Init(NULL, _vq, NULL); }
   VectorL2MortarIntegrator(MatrixCoefficient *_mq) { Init(NULL, NULL, _mq); }

   void AssembleElementMatrix(const FiniteElement &trial,
                              const IntegrationRule &trial_ir,
                              ElementTransformation &trial_Trans,
                              const FiniteElement &test,
                              const IntegrationRule &test_ir,
                              ElementTransformation &test_Trans,
                              DenseMatrix &elemmat) override;

   inline bool is_vector_fe() const override { return true; }

private:
   Coefficient *Q;
   VectorCoefficient *VQ;
   MatrixCoefficient *MQ;

   void Init(Coefficient *q, VectorCoefficient *vq, MatrixCoefficient *mq)
   {
      Q = q;
      VQ = vq;
      MQ = mq;
   }
};

} // namespace mfem

#endif // MFEML2_MORTAR_INTEGRATOR_HPP