#pragma once

#include "CSolver.hpp"
#include "../variables/CDiscAdjFEABoundVariable.hpp"

class CDiscAdjFEASolver final : public CSolver {
private:
  unsigned short KindDirect_Solver = 0;
  CSolver *direct_solver = nullptr;
  su2double *Sens_E = nullptr,          
  *Sens_Nu = nullptr,                   
  *Sens_nL = nullptr,                   
  **CSensitivity = nullptr;             

  su2double *Solution_Vel = nullptr,    
  *Solution_Accel = nullptr;            

  su2double ObjFunc_Value = 0.0;        
  su2double *normalLoads = nullptr;     

  unsigned short nMPROP = 0;            

  su2double *E_i = nullptr,               
            *Nu_i = nullptr,              
            *Rho_i = nullptr,             
            *Rho_DL_i = nullptr;          
  int       *AD_Idx_E_i = nullptr,        
            *AD_Idx_Nu_i = nullptr,       
            *AD_Idx_Rho_i = nullptr,      
            *AD_Idx_Rho_DL_i = nullptr;   

  su2double *Local_Sens_E = nullptr,        
            *Global_Sens_E = nullptr,       
            *Total_Sens_E = nullptr;        
  su2double *Local_Sens_Nu = nullptr,       
            *Global_Sens_Nu = nullptr,      
            *Total_Sens_Nu = nullptr;       
  su2double *Local_Sens_Rho = nullptr,      
            *Global_Sens_Rho = nullptr,     
            *Total_Sens_Rho = nullptr;      
  su2double *Local_Sens_Rho_DL = nullptr,   
            *Global_Sens_Rho_DL = nullptr,  
            *Total_Sens_Rho_DL = nullptr;   

  bool de_effects = false;                  
  unsigned short nEField = 0;               
  su2double *EField = nullptr;              
  int       *AD_Idx_EField = nullptr;       
  su2double *Local_Sens_EField = nullptr,   
            *Global_Sens_EField = nullptr,  
            *Total_Sens_EField = nullptr;   

  bool fea_dv = false;                  
  unsigned short nDV = 0;               
  su2double *DV_Val = nullptr;          
  int       *AD_Idx_DV_Val = nullptr;   
  su2double *Local_Sens_DV = nullptr,   
            *Global_Sens_DV = nullptr,  
            *Total_Sens_DV = nullptr;   

  CDiscAdjFEABoundVariable* nodes = nullptr;  

  inline CVariable* GetBaseClassPointerToNodes() override { return nodes; }

public:

  CDiscAdjFEASolver(void);

  CDiscAdjFEASolver(CGeometry *geometry, CConfig *config);

  CDiscAdjFEASolver(CGeometry *geometry, CConfig *config, CSolver* solver, unsigned short Kind_Solver, unsigned short iMesh);

  ~CDiscAdjFEASolver(void) override;

  void RegisterSolution(CGeometry *geometry, CConfig *config) override;

  void RegisterOutput(CGeometry *geometry, CConfig *config) override;

  void SetAdjoint_Output(CGeometry *geometry, CConfig *config) override;

  void ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config) override;

  void RegisterObj_Func(CConfig *config) override;

  void SetSurface_Sensitivity(CGeometry *geometry, CConfig* config) override;

  void SetSensitivity(CGeometry *geometry, CSolver **solver, CConfig *config) override;

  void SetAdj_ObjFunc(CGeometry *geometry, CConfig* config) override;

  inline su2double GetTotal_Sens_E(unsigned short iVal) const override { return Total_Sens_E[iVal]; }

  inline su2double GetTotal_Sens_Nu(unsigned short iVal) const override { return Total_Sens_Nu[iVal]; }

  inline su2double GetTotal_Sens_Rho(unsigned short iVal) const override { return Total_Sens_Rho[iVal]; }

  inline su2double GetTotal_Sens_Rho_DL(unsigned short iVal) const override { return Total_Sens_Rho_DL[iVal]; }

  inline su2double GetTotal_Sens_EField(unsigned short iEField) const override { return Total_Sens_EField[iEField]; }

  inline su2double GetTotal_Sens_DVFEA(unsigned short iDVFEA) const override { return Total_Sens_DV[iDVFEA]; }

  inline su2double GetGlobal_Sens_E(unsigned short iVal) const override { return Global_Sens_E[iVal]; }

  inline su2double GetGlobal_Sens_Nu(unsigned short iVal) const override { return Global_Sens_Nu[iVal]; }

  inline su2double GetGlobal_Sens_EField(unsigned short iEField) const override { return Global_Sens_EField[iEField]; }

  inline su2double GetGlobal_Sens_DVFEA(unsigned short iDVFEA) const override { return Global_Sens_DV[iDVFEA]; }

  inline su2double GetGlobal_Sens_Rho(unsigned short iVal) const override { return Global_Sens_Rho[iVal]; }

  inline su2double GetGlobal_Sens_Rho_DL(unsigned short iVal) const override { return Global_Sens_Rho_DL[iVal]; }

  inline su2double GetVal_Young(unsigned short iVal) const override { return E_i[iVal]; }

  inline su2double GetVal_Poisson(unsigned short iVal) const override { return Nu_i[iVal]; }

  inline su2double GetVal_Rho(unsigned short iVal) const override { return Rho_i[iVal]; }

  inline su2double GetVal_Rho_DL(unsigned short iVal) const override { return Rho_DL_i[iVal]; }

  inline unsigned short GetnEField(void) const override { return nEField; }

  void ReadDV(CConfig *config) override;

  inline unsigned short GetnDVFEA(void) const override { return nDV; }

  inline su2double GetVal_EField(unsigned short iVal) const override { return EField[iVal]; }

  inline su2double GetVal_DVFEA(unsigned short iVal) const override { return DV_Val[iVal]; }

  void SetRecording(CGeometry *geometry, CConfig *config) override ;

  void RegisterVariables(CGeometry *geometry,
                         CConfig *config,
                         bool reset = false) override;

  void ExtractAdjoint_Variables(CGeometry *geometry, CConfig *config) override;

  void Preprocessing(CGeometry *geometry,
                    CSolver **solver_container,
                    CConfig *config,
                    unsigned short iMesh,
                    unsigned short iRKStep,
                    unsigned short RunTime_EqSystem,
                    bool Output) override;

  void LoadRestart(CGeometry **geometry,
                   CSolver ***solver,
                   CConfig *config,
                   int val_iter,
                   bool val_update_geo) override;

};
/*!
 * \file CDiscAdjFEASolver.hpp
 * \brief Headers of the CDiscAdjFEASolver class
 * \author R. Sanchez
 * \version 7.0.6 "Blackbird"
 *
 * SU2 Project Website: https://su2code.github.io
 *
 * The SU2 Project is maintained by the SU2 Foundation
 * (http://su2foundation.org)
 *
 * Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md)
 *
 * SU2 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.
 *
 * SU2 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 SU2. If not, see <http://www.gnu.org/licenses/>.
 */

#pragma once

#include "CSolver.hpp"
#include "../variables/CDiscAdjFEABoundVariable.hpp"

/*!
 * \class CDiscAdjFEASolver
 * \brief Main class for defining the discrete adjoint solver for FE structural problems.
 * \ingroup Discrete_Adjoint
 * \author R. Sanchez
 */
class CDiscAdjFEASolver final : public CSolver {
private:
  unsigned short KindDirect_Solver = 0;
  CSolver *direct_solver = nullptr;
  su2double *Sens_E = nullptr,          /*!< \brief Young modulus sensitivity coefficient for each boundary. */
  *Sens_Nu = nullptr,                   /*!< \brief Poisson's ratio sensitivity coefficient for each boundary. */
  *Sens_nL = nullptr,                   /*!< \brief Normal pressure sensitivity coefficient for each boundary. */
  **CSensitivity = nullptr;             /*!< \brief Shape sensitivity coefficient for each boundary and vertex. */

  su2double *Solution_Vel = nullptr,    /*!< \brief Velocity componenent of the solution. */
  *Solution_Accel = nullptr;            /*!< \brief Acceleration componenent of the solution. */

  su2double ObjFunc_Value = 0.0;        /*!< \brief Value of the objective function. */
  su2double *normalLoads = nullptr;     /*!< \brief Values of the normal loads for each marker iMarker_nL. */

  unsigned short nMPROP = 0;            /*!< \brief Number of material properties */

  su2double *E_i = nullptr,               /*!< \brief Values of the Young's Modulus. */
            *Nu_i = nullptr,              /*!< \brief Values of the Poisson's ratio. */
            *Rho_i = nullptr,             /*!< \brief Values of the density (for inertial effects). */
            *Rho_DL_i = nullptr;          /*!< \brief Values of the density (for volume loading). */
  int       *AD_Idx_E_i = nullptr,        /*!< \brief Derivative index of the Young's Modulus. */
            *AD_Idx_Nu_i = nullptr,       /*!< \brief Derivative index of the Poisson's ratio. */
            *AD_Idx_Rho_i = nullptr,      /*!< \brief Derivative index of the density (for inertial effects). */
            *AD_Idx_Rho_DL_i = nullptr;   /*!< \brief Derivative index of the density (for volume loading). */

  su2double *Local_Sens_E = nullptr,        /*!< \brief Local sensitivity of the Young's modulus. */
            *Global_Sens_E = nullptr,       /*!< \brief Global sensitivity of the Young's modulus. */
            *Total_Sens_E = nullptr;        /*!< \brief Total sensitivity of the Young's modulus (time domain). */
  su2double *Local_Sens_Nu = nullptr,       /*!< \brief Local sensitivity of the Poisson ratio. */
            *Global_Sens_Nu = nullptr,      /*!< \brief Global sensitivity of the Poisson ratio. */
            *Total_Sens_Nu = nullptr;       /*!< \brief Total sensitivity of the Poisson ratio (time domain). */
  su2double *Local_Sens_Rho = nullptr,      /*!< \brief Local sensitivity of the density. */
            *Global_Sens_Rho = nullptr,     /*!< \brief Global sensitivity of the density. */
            *Total_Sens_Rho = nullptr;      /*!< \brief Total sensitivity of the density (time domain). */
  su2double *Local_Sens_Rho_DL = nullptr,   /*!< \brief Local sensitivity of the volume load. */
            *Global_Sens_Rho_DL = nullptr,  /*!< \brief Global sensitivity of the volume load. */
            *Total_Sens_Rho_DL = nullptr;   /*!< \brief Total sensitivity of the volume load (time domain). */

  bool de_effects = false;                  /*!< \brief Determines if DE effects are considered. */
  unsigned short nEField = 0;               /*!< \brief Number of electric field areas in the problem. */
  su2double *EField = nullptr;              /*!< \brief Array that stores the electric field as design variables. */
  int       *AD_Idx_EField = nullptr;       /*!< \brief Derivative index of the electric field as design variables. */
  su2double *Local_Sens_EField = nullptr,   /*!< \brief Local sensitivity of the Electric Field. */
            *Global_Sens_EField = nullptr,  /*!< \brief Global sensitivity of the Electric Field. */
            *Total_Sens_EField = nullptr;   /*!< \brief Total sensitivity of the Electric Field (time domain). */

  bool fea_dv = false;                  /*!< \brief Determines if the design variable we study is a FEA parameter. */
  unsigned short nDV = 0;               /*!< \brief Number of design variables in the problem. */
  su2double *DV_Val = nullptr;          /*!< \brief Values of the design variables. */
  int       *AD_Idx_DV_Val = nullptr;   /*!< \brief Derivative index of the design variables. */
  su2double *Local_Sens_DV = nullptr,   /*!< \brief Local sensitivity of the design variables. */
            *Global_Sens_DV = nullptr,  /*!< \brief Global sensitivity of the design variables. */
            *Total_Sens_DV = nullptr;   /*!< \brief Total sensitivity of the design variables (time domain). */

  CDiscAdjFEABoundVariable* nodes = nullptr;  /*!< \brief The highest level in the variable hierarchy this solver can safely use. */

  /*!
   * \brief Return nodes to allow CSolver::base_nodes to be set.
   */
  inline CVariable* GetBaseClassPointerToNodes() override { return nodes; }

public:

  /*!
   * \brief Constructor of the class.
   */
  CDiscAdjFEASolver(void);

  /*!
   * \overload
   * \param[in] geometry - Geometrical definition of the problem.
   * \param[in] config - Definition of the particular problem.
   * \param[in] iMesh - Index of the mesh in multigrid computations.
   */
  CDiscAdjFEASolver(CGeometry *geometry, CConfig *config);

  /*!
   * \overload
   * \param[in] geometry - Geometrical definition of the problem.
   * \param[in] config - Definition of the particular problem.
   * \param[in] solver - Initialize the discrete adjoint solver with the corresponding direct solver.
   * \param[in] Kind_Solver - The kind of direct solver.
   */
  CDiscAdjFEASolver(CGeometry *geometry, CConfig *config, CSolver* solver, unsigned short Kind_Solver, unsigned short iMesh);

  /*!
   * \brief Destructor of the class.
   */
  ~CDiscAdjFEASolver(void) override;

  /*!
   * \brief Performs the preprocessing of the adjoint AD-based solver.
   *        Registers all necessary variables on the tape. Called while tape is active.
   * \param[in] geometry_container - The geometry container holding all grid levels.
   * \param[in] config_container - The particular config.
   */
  void RegisterSolution(CGeometry *geometry, CConfig *config) override;

  /*!
   * \brief Performs the preprocessing of the adjoint AD-based solver.
   *        Registers all necessary variables that are output variables on the tape.
   *        Called while tape is active.
   * \param[in] geometry_container - The geometry container holding all grid levels.
   * \param[in] config_container - The particular config.
   */
  void RegisterOutput(CGeometry *geometry, CConfig *config) override;

  /*!
   * \brief Sets the adjoint values of the output of the flow (+turb.) iteration
   *         before evaluation of the tape.
   * \param[in] geometry - The geometrical definition of the problem.
   * \param[in] config - The particular config.
   */
  void SetAdjoint_Output(CGeometry *geometry, CConfig *config) override;

  /*!
   * \brief Sets the adjoint values of the input variables of the flow (+turb.) iteration
   *        after tape has been evaluated.
   * \param[in] geometry - The geometrical definition of the problem.
   * \param[in] config - The particular config.
   */
  void ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config) override;

  /*!
   * \brief Register the objective function as output.
   * \param[in] geometry - The geometrical definition of the problem.
   */
  void RegisterObj_Func(CConfig *config) override;

  /*!
   * \brief Set the surface sensitivity.
   * \param[in] geometry - Geometrical definition of the problem.
   * \param[in] config - Definition of the particular problem.
   */
  void SetSurface_Sensitivity(CGeometry *geometry, CConfig* config) override;

  /*!
   * \brief Extract and set the geometrical sensitivity.
   * \param[in] geometry - Geometrical definition of the problem.
   * \param[in] solver - The solver container holding all terms of the solution.
   * \param[in] config - Definition of the particular problem.
   */
  void SetSensitivity(CGeometry *geometry, CSolver **solver, CConfig *config) override;

  /*!
   * \brief Set the objective function.
   * \param[in] geometry - Geometrical definition of the problem.
   * \param[in] config - Definition of the particular problem.
   */
  void SetAdj_ObjFunc(CGeometry *geometry, CConfig* config) override;

  /*!
   * \brief Provide the total Young's modulus sensitivity
   * \return Value of the total Young's modulus sensitivity
   *         (inviscid + viscous contribution).
   */
  inline su2double GetTotal_Sens_E(unsigned short iVal) const override { return Total_Sens_E[iVal]; }

  /*!
   * \brief Set the total Poisson's ratio sensitivity.
   * \return Value of the Poisson's ratio sensitivity
   */
  inline su2double GetTotal_Sens_Nu(unsigned short iVal) const override { return Total_Sens_Nu[iVal]; }

  /*!
   * \brief Get the total sensitivity for the structural density
   * \return Value of the structural density sensitivity
   */
  inline su2double GetTotal_Sens_Rho(unsigned short iVal) const override { return Total_Sens_Rho[iVal]; }

  /*!
   * \brief Get the total sensitivity for the structural weight
   * \return Value of the structural weight sensitivity
   */
  inline su2double GetTotal_Sens_Rho_DL(unsigned short iVal) const override { return Total_Sens_Rho_DL[iVal]; }

  /*!
   * \brief A virtual member.
   * \return Value of the sensitivity coefficient for the Electric Field in the region iEField (time averaged)
   */
  inline su2double GetTotal_Sens_EField(unsigned short iEField) const override { return Total_Sens_EField[iEField]; }

  /*!
   * \brief A virtual member.
   * \return Value of the total sensitivity coefficient for the FEA DV in the region iDVFEA (time averaged)
   */
  inline su2double GetTotal_Sens_DVFEA(unsigned short iDVFEA) const override { return Total_Sens_DV[iDVFEA]; }

  /*!
   * \brief A virtual member.
   * \return Value of the sensitivity coefficient for the Young Modulus E
   */
  inline su2double GetGlobal_Sens_E(unsigned short iVal) const override { return Global_Sens_E[iVal]; }

  /*!
   * \brief A virtual member.
   * \return Value of the Mach sensitivity for the Poisson's ratio Nu
   */
  inline su2double GetGlobal_Sens_Nu(unsigned short iVal) const override { return Global_Sens_Nu[iVal]; }

  /*!
   * \brief A virtual member.
   * \return Value of the sensitivity coefficient for the Electric Field in the region iEField
   */
  inline su2double GetGlobal_Sens_EField(unsigned short iEField) const override { return Global_Sens_EField[iEField]; }

  /*!
   * \brief A virtual member.
   * \return Value of the sensitivity coefficient for the FEA DV in the region iDVFEA
   */
  inline su2double GetGlobal_Sens_DVFEA(unsigned short iDVFEA) const override { return Global_Sens_DV[iDVFEA]; }

  /*!
   * \brief Get the total sensitivity for the structural density
   * \return Value of the structural density sensitivity
   */
  inline su2double GetGlobal_Sens_Rho(unsigned short iVal) const override { return Global_Sens_Rho[iVal]; }

  /*!
   * \brief Get the total sensitivity for the structural weight
   * \return Value of the structural weight sensitivity
   */
  inline su2double GetGlobal_Sens_Rho_DL(unsigned short iVal) const override { return Global_Sens_Rho_DL[iVal]; }

  /*!
   * \brief Get the value of the Young modulus from the adjoint solver
   * \return Value of the Young modulus from the adjoint solver
   */
  inline su2double GetVal_Young(unsigned short iVal) const override { return E_i[iVal]; }

  /*!
   * \brief Get the value of the Poisson's ratio from the adjoint solver
   * \return Value of the Poisson's ratio from the adjoint solver
   */
  inline su2double GetVal_Poisson(unsigned short iVal) const override { return Nu_i[iVal]; }

  /*!
   * \brief Get the value of the density from the adjoint solver, for inertial effects
   * \return Value of the density from the adjoint solver
   */
  inline su2double GetVal_Rho(unsigned short iVal) const override { return Rho_i[iVal]; }

  /*!
   * \brief Get the value of the density from the adjoint solver, for dead loads
   * \return Value of the density for dead loads, from the adjoint solver
   */
  inline su2double GetVal_Rho_DL(unsigned short iVal) const override { return Rho_DL_i[iVal]; }

  /*!
   * \brief Get the number of variables for the Electric Field from the adjoint solver
   * \return Number of electric field variables from the adjoint solver
   */
  inline unsigned short GetnEField(void) const override { return nEField; }

  /*!
   * \brief Read the design variables for the adjoint solver
   */
  void ReadDV(CConfig *config) override;

  /*!
   * \brief Get the number of design variables from the adjoint solver,
   * \return Number of design variables from the adjoint solver
   */
  inline unsigned short GetnDVFEA(void) const override { return nDV; }

  /*!
   * \brief Get the value of the Electric Field from the adjoint solver
   * \return Pointer to the values of the Electric Field
   */
  inline su2double GetVal_EField(unsigned short iVal) const override { return EField[iVal]; }

  /*!
   * \brief Get the value of the design variables from the adjoint solver
   * \return Pointer to the values of the design variables
   */
  inline su2double GetVal_DVFEA(unsigned short iVal) const override { return DV_Val[iVal]; }

  /*!
   * \brief Prepare the solver for a new recording.
   * \param[in] kind_recording - Kind of AD recording.
   */
  void SetRecording(CGeometry *geometry, CConfig *config) override ;

  /*!
   * \brief A virtual member.
   * \param[in] geometry - Geometrical definition of the problem.
   * \param[in] config - Definition of the particular problem.
   * \param[in] reset - If true reset variables to their initial values.
   */
  void RegisterVariables(CGeometry *geometry,
                         CConfig *config,
                         bool reset = false) override;

  /*!
   * \brief A virtual member.
   * \param[in] geometry - Geometrical definition of the problem.
   * \param[in] config - Definition of the particular problem.
   */
  void ExtractAdjoint_Variables(CGeometry *geometry, CConfig *config) override;

  /*!
   * \brief Update the dual-time derivatives.
   * \param[in] geometry - Geometrical definition of the problem.
   * \param[in] solver_container - Container vector with all the solutions.
   * \param[in] config - Definition of the particular problem.
   * \param[in] iMesh - Index of the mesh in multigrid computations.
   * \param[in] iRKStep - Current step of the Runge-Kutta iteration.
   * \param[in] RunTime_EqSystem - System of equations which is going to be solved.
   * \param[in] Output - boolean to determine whether to print output.
   */
  void Preprocessing(CGeometry *geometry,
                    CSolver **solver_container,
                    CConfig *config,
                    unsigned short iMesh,
                    unsigned short iRKStep,
                    unsigned short RunTime_EqSystem,
                    bool Output) override;

  /*!
   * \brief Load a solution from a restart file.
   * \param[in] geometry - Geometrical definition of the problem.
   * \param[in] solver - Container vector with all of the solvers.
   * \param[in] config - Definition of the particular problem.
   * \param[in] val_iter - Current external iteration number.
   * \param[in] val_update_geo - Flag for updating coords and grid velocity.
   */
  void LoadRestart(CGeometry **geometry,
                   CSolver ***solver,
                   CConfig *config,
                   int val_iter,
                   bool val_update_geo) override;

};
