// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL  : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program 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 General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.

//#include "HyperDN.h"

     // Constructeur par defaut
template <class TensHH,class TensBB,class TensBH,class TensHB,
          class Tens_nHH,class Tens_nBB>
HyperDN<TensHH,TensBB,TensBH,TensHB,Tens_nHH,Tens_nBB>::HyperDN () :
  Loi_comp_abstraite()
	 {}
	// Constructeur utile si l'identificateur du nom de la loi
	// de comportement et la dimension  sont connus
template <class TensHH,class TensBB,class TensBH,class TensHB,
          class Tens_nHH,class Tens_nBB>
HyperDN<TensHH,TensBB,TensBH,TensHB,Tens_nHH,Tens_nBB>
             ::HyperDN (Enum_comp id_compor,Enum_categorie_loi_comp categorie_comp,int dimension) :
	 Loi_comp_abstraite(id_compor,categorie_comp,dimension)
	  {} 
	// Constructeur utile si l'identificateur du nom de la loi
	// de comportement et la dimension  sont connus			
template <class TensHH,class TensBB,class TensBH,class TensHB,
          class Tens_nHH,class Tens_nBB>
HyperDN<TensHH,TensBB,TensBH,TensHB,Tens_nHH,Tens_nBB>
             ::HyperDN (char* nom,Enum_categorie_loi_comp categorie_comp,int dimension) :
      Loi_comp_abstraite(nom,categorie_comp,dimension)
	  {} 	// DESTRUCTEUR :
template <class TensHH,class TensBB,class TensBH,class TensHB,
          class Tens_nHH,class Tens_nBB>
HyperDN<TensHH,TensBB,TensBH,TensHB,Tens_nHH,Tens_nBB>
             ::~HyperDN ()
	  {} 
	// constructeur de copie
template <class TensHH,class TensBB,class TensBH,class TensHB,
          class Tens_nHH,class Tens_nBB>
HyperDN<TensHH,TensBB,TensBH,TensHB,Tens_nHH,Tens_nBB>
             ::HyperDN (const HyperDN & a) :
	 Loi_comp_abstraite (a)
	  {}
	  
  
 // ========== codage des METHODES VIRTUELLES  protegees:================
        // calcul des contraintes a tdt

template <class TensHH,class TensBB,class TensBH,class TensHB,
          class Tens_nHH,class Tens_nBB>
void HyperDN<TensHH,TensBB,TensBH,TensHB,Tens_nHH,Tens_nBB>
              ::Calcul_SigmaHH 
             (TenseurHH& ,TenseurBB& ,DdlElement & tab_ddl,
             TenseurBB & ,TenseurHH & ,BaseB& ,BaseH& ,TenseurBB & epsBB_,TenseurBB & ,
             TenseurBB & gijBB_,TenseurHH & gijHH_,Tableau <TenseurBB *>& d_gijBB_,
             double& jacobien_0,double& jacobien_,TenseurHH & sigHH
		  	,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double&  module_cisaillement
		  	,const Met_abstraite::Expli_t_tdt&  )
 {   
   #ifdef MISE_AU_POINT	 	 
	 if (tab_ddl.NbDdl() != d_gijBB_.Taille())
	    { cout << "\nErreur : le nb de ddl est != de la taille de d_gijBB_  !\n";
		  cout << " HyperDN::Calcul_SigmaHH\n";
		  Sortie(1);
		};		
    #endif
    TensBH epsBH; // la déformation en mixte
    TensBH * IdGBH; // l'identitée
    double  alpha_0,alpha_1,alpha_2; // coefficients pour le calcul de la loi de comportement
    // Calcul des 3 invariants , de epsBH, 
    // retour de IdGBH qui pointe sur le bon tenseur
    double Ieps,V,bIIb; 
    IdGBH = Invariants (epsBB_,gijBB_,gijHH_,jacobien_0,jacobien_,Ieps,V,bIIb,epsBH);
    // calcul du potentiel et de ses dérivées premières
    double E,EV,EbIIb,EIeps; // potentiel et dérivées
    Potentiel(jacobien_0,Ieps,V,bIIb,E,EV,EbIIb,EIeps); 
    // calcul des alphas
    Alpha(E,EV,EbIIb,EIeps,jacobien_0,Ieps,V,bIIb,alpha_0,alpha_1,alpha_2);
    // calcul des contraintes finales     
    TensBH  sigBH = alpha_0 * (*IdGBH) + alpha_1 * epsBH 
                   + alpha_2 * (epsBH * epsBH);
    sigHH = gijHH_ * sigBH;  // en deux fois contravariants   
 }
        
 // calcul des contraintes et de ses variations  a t+dt
template <class TensHH,class TensBB,class TensBH,class TensHB,
          class Tens_nHH,class Tens_nBB>
 void HyperDN<TensHH,TensBB,TensBH,TensHB,Tens_nHH,Tens_nBB>
          ::Calcul_DsigmaHH_tdt 
           (TenseurHH& ,TenseurBB& ,DdlElement & tab_ddl
           ,BaseB& ,TenseurBB & ,TenseurHH & ,
            BaseB& ,Tableau <BaseB> & ,BaseH& ,Tableau <BaseH> & ,
            TenseurBB & epsBB_tdt,Tableau <TenseurBB *>& d_epsBB,
            TenseurBB & ,TenseurBB & gijBB_tdt,TenseurHH & gijHH_tdt,
		    Tableau <TenseurBB *>& d_gijBB_tdt,
		  	Tableau <TenseurHH *>& d_gijHH_tdt,double& jacobien_0,double& jacobien_tdt,
		  	Vecteur& d_jacobien_tdt,TenseurHH& sigHH,Tableau <TenseurHH *>& d_sigHH
		  	,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double&  module_cisaillement
		  	,const Met_abstraite::Impli& )

 {   
   #ifdef MISE_AU_POINT	 	 
	 if (tab_ddl.NbDdl() != d_gijBB_tdt.Taille())
	    { cout << "\nErreur : le nb de ddl est != de la taille de d_gijBB_t  !\n";
		  cout << " HyperDN::Calcul_DsigmaHH_tdt\n";
		  Sortie(1);
		};		
    #endif
    int nbddl = tab_ddl.NbDdl();
    TensBH epsBH; // la déformation en mixte
    Tableau<TensBH>  depsBH(nbddl);
    TensBH * IdGBH; // l'identitée
    // def des alphas
    double  alpha_0,alpha_1,alpha_2; 
    Tableau<double> dalpha_0(nbddl),dalpha_1(nbddl),dalpha_2(nbddl); 
    // calcul des trois invariants et de leurs variations, de epsBH, 
    // et de sa variation, puis retour de IdGBH qui pointe sur le bon tenseur
    double Ieps,V,bIIb;
    Tableau<double> dIeps(nbddl),dV(nbddl),dbIIb(nbddl); 
    IdGBH = Invariants_et_var (epsBB_tdt,gijBB_tdt,d_gijBB_tdt,
           gijHH_tdt,d_gijHH_tdt,jacobien_0,jacobien_tdt,d_jacobien_tdt,
           Ieps,dIeps,V,dV,bIIb,dbIIb,epsBH,depsBH);
    // calcul  du potentiel et de ses dérivées premières et secondes
    double E,EV,EbIIb,EIeps,EVV,EbIIb2,EIeps2,EVbIIb,EVIeps,EbIIbIeps; // potentiel et variations
    Potentiel_et_var(jacobien_0,Ieps,V,bIIb,E,EV,EbIIb,EIeps,EVV,EbIIb2,EIeps2,
                  EVbIIb,EVIeps,EbIIbIeps );
    // calcul des coefficients alpha  et de leurs variations par rapport aux degrés de libertés
    Tableau<double> dE(nbddl),dEV(nbddl),dEbIIb(nbddl),dEVV(nbddl),dEQQ(nbddl),dEVQ(nbddl);
    Alpha_var ( E, EV, EbIIb, EIeps,EVV, EbIIb2, EIeps2, EVbIIb, EVIeps, EbIIbIeps,
              jacobien_0,Ieps,dIeps,V,dV,bIIb,dbIIb,alpha_0,dalpha_0,alpha_1,dalpha_1,alpha_2,dalpha_2);
    // calcul des contraintes finales     
    TensBH epsepsBH = epsBH * epsBH; // epsBH:epsBH
    TensBH  sigBH = alpha_0 * (*IdGBH) + alpha_1 * epsBH 
                   + alpha_2 * epsepsBH;
    sigHH = gijHH_tdt * sigBH;  // en deux fois contravariants  
     
    // cas le la variation du tenseur des contraintes
        
    for (int i = 1; i<= nbddl; i++)
     { const TensHH & dgijHH = (*(d_gijHH_tdt(i))) ;        // pour simplifier l'ecriture
       TensHH& dsigHH = *((TensHH*) (d_sigHH(i)));      //
       
       TensBB depsBB = 0.5 * (*(d_gijBB_tdt(i)));
       TensBH depsBH =  epsBB_tdt * dgijHH  + depsBB * gijHH_tdt ;
       TensBH  dsigBH = dalpha_0(i) * (*IdGBH) + dalpha_1(i) * epsBH 
                   + dalpha_2(i) * epsepsBH + alpha_1 * depsBH + alpha_2 * (depsBH*epsBH + epsBH*depsBH);
       dsigHH = dgijHH * sigBH + gijHH_tdt * dsigBH;        
      }          
  }

// -------------------------------------------------------------------------------------------------------------
    // METHODES PROTEGEES :
// -------------------------------------------------------------------------------------------------------------
    
 // calcul des coefficients alpha 
template <class TensHH,class TensBB,class TensBH,class TensHB,
          class Tens_nHH,class Tens_nBB>
 inline void HyperDN<TensHH,TensBB,TensBH,TensHB,Tens_nHH,Tens_nBB>
           ::Alpha
 (double& ,double& EV,double& EbIIb,double& EIeps,
        double& jacobien_0,double & Ieps,double & V,double& ,
        double & alpha_0,double & alpha_1,double & alpha_2)
  { 
    double unsurrgOV = 1./(jacobien_0*V);
    alpha_0 = unsurrgOV * ( EIeps + V * EV - Ieps/3. * EbIIb);
    alpha_1 = unsurrgOV * ( -2. * EIeps  + (2 * Ieps/3. + 1.) * EbIIb);
    alpha_2 = unsurrgOV * (  - 2. * EbIIb);
   }
     
 // calcul des coefficients alpha  et de leurs variations par rapport aux degrés de libertés
template <class TensHH,class TensBB,class TensBH,class TensHB,
          class Tens_nHH,class Tens_nBB>
 inline void HyperDN<TensHH,TensBB,TensBH,TensHB,Tens_nHH,Tens_nBB>
           :: Alpha_var (
        double& ,double& EV,double& EbIIb,double& EIeps,
        double& EVV,double& EbIIb2,double& EIeps2,
        double& EVbIIb,double& EVIeps,double& EbIIbIeps,
        double& jacobien_0,double&  Ieps,Tableau<double> & dIeps,double&  V,Tableau<double> & dV,
        double& ,Tableau<double> & dbIIb,
        double&  alpha_0,Tableau<double> & dalpha_0,
        double&  alpha_1,Tableau<double> & dalpha_1,
        double&  alpha_2,Tableau<double> & dalpha_2)
  { // calcul des coefficients alpha
    double  unsurrgOV = 1./(jacobien_0*V);
    double al1 = ( EIeps + V * EV - Ieps/3. * EbIIb); // pour optimiser les calculs de variations
    alpha_0 = unsurrgOV * al1;
    double al2 = ( -2. * EIeps  + (2 * Ieps/3. + 1.) * EbIIb); // pour optimiser les calculs de variations
    alpha_1 = unsurrgOV * al2;
    double al3 = (  - 2. * EbIIb); // pour optimiser les calculs de variations
    alpha_2 = unsurrgOV * al3;
    // cas le la variation des coefficients alpha
    int nbddl = dIeps.Taille();
    for (int i=1; i<= nbddl; i++)
      {  // variation de  unsurrgOV
         double d_unsurrgOV = -  unsurrgOV* unsurrgOV * (jacobien_0 * dV(i));
         // variations des alphas
         dalpha_0(i) = d_unsurrgOV * al1 + unsurrgOV *
          ( dIeps(i) * EIeps2 + dbIIb(i) * EbIIbIeps + dV(i) * EVIeps);
         dalpha_1(i) = d_unsurrgOV * al2 + unsurrgOV *
          ( dIeps(i) * EbIIbIeps + dbIIb(i) * EbIIb2 + dV(i) * EVbIIb);          
         dalpha_2(i) = d_unsurrgOV * al3 + unsurrgOV *
          ( dIeps(i) * EVIeps + dbIIb(i) * EVbIIb + dV(i) * EVV);
        };   
   }