// 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 "Debug.h"
# include "Deformation.h"
#include "ConstMath.h"
#include "ParaGlob.h"
#include "MathUtil.h"
#include "Tenseur3.h"
#include "Tenseur2.h"
#include "Tenseur1.h"
# include "TypeConsTens.h"

// -----------------------------------------------------------------
//      cas de la déformation d'Almansi
// -----------------------------------------------------------------

  	  // cas implicite : tous les parametres sont de resultats	
const Met_abstraite::Impli& Deformation::Cal_implicit_Almansi (bool 
           ,TenseurBB & epsBB_tdt,Tableau <TenseurBB *> & d_epsBB_tdt
           ,TenseurBB& DepsBB,TenseurBB& delta_epsBB,bool premier_calcul
           ,const Met_abstraite::Impli& ex)
  { 
     const VariablesTemps& vartemps = ParaGlob::Variables_de_temps();
     double deltat=vartemps.IncreTempsCourant();
// modif du 3 nov 2014 -> introduction de la sauvegarde de D : D_BB_t et D_BB_tdt
//     // dans le cas où le pas de temps est trop petit on met un message d'erreur
//     if (Dabs(deltat) < ConstMath::trespetit)
//     	{ cout << "\n erreur le pas de temps " << deltat
//     	       << " est trop faible pour un calcul correcte de la vitesse moyenne";
//     	  if (ParaGlob::NiveauImpression() > 3)
//     	       cout << "\n Deformation::Cal_implicit_Almansi ( .... ";
//     	  Sortie(1);
//     	};
     // dans le cas où ce n'est pas le premier calcul on récupère les datas qui ne sont pas recalculée
     if (!premier_calcul)
      { Deformation::Stmet& a_0 = saveDefResul->meti_00;  // pour simplifier
        Deformation::Stmet& a_t = saveDefResul->meti_t;  // pour simplifier
        metrique->Recup_grandeur_0_t(*a_0.giB_,*a_0.giH_,*a_t.giB_,*a_t.giH_
                                     ,*a_0.gijBB_,*a_0.gijHH_,*a_t.gijBB_,*a_t.gijHH_
                                     ,*a_t.gradVmoyBB_,*a_t.jacobien_,*a_0.jacobien_);
      };
////--- debug
//cout << "\n Deformation::Cal_implicit_Almansi ";
//cout << "\n gijBB_t= ";ex.gijBB_t->Ecriture(cout);
//cout << "\n stocké meti_t.gijBB_ ";saveDefResul->meti_t.gijBB_->Ecriture(cout);
////--- fin debug
   
     // -- calcul de la déformation d'Almansi et de l'incrément de déformation
     //     epsBB_tdt et delta_epsBB est dimensionne auparavant             
     epsBB_tdt = 0.5 * (*(ex.gijBB_tdt) - *(ex.gijBB_0));
     delta_epsBB = 0.5 * (*(ex.gijBB_tdt) - *(ex.gijBB_t));
     // -- calcul de la vitesse de déformation: a priori on utilise 
     // la variation du tenseur déformation sur delta t ce qui correspond à l'utilisation du gradient
     // de vitesse de déformation moyen calculé à t+(delta t)/2
     if (Dabs(deltat) >= ConstMath::trespetit)
        { DepsBB = delta_epsBB/deltat; }
// modif du 3 nov 2014 -> introduction de la sauvegarde de D : D_BB_t et D_BB_tdt
//     else // dans le cas où l'incrément de temps est nul pour l'instant on met la vitesse de déformation à 0
//        { DepsBB.Inita(0.);}
     else
        // on utilise la vitesse précédente pour la vitesse de déformation
        { DepsBB = (*saveDefResul->D_BB_t);};
   
//----------------- debug
//TenseurBB truc =deltat*DepsBB;
//cout << "\n ------------- du cote deformation ---------------";
//cout << "\n epsBB="; epsBB_tdt.Ecriture(cout);
//cout << "\n delta_epsBB=";delta_epsBB.Ecriture(cout);
//
//cout << "\n DepsBB*deltat="<<truc<<" ";
//----------- fin debug        
	 
//	    //on utilise le gradient de vitesse moyen -> pas bon car le gradient de vitesse moyen sera également null donc calcul inutile
//	      // pour éviter une division par 0
//	   { DepsBB = 0.5 * ((*ex.gradVmoyBB_tdt) + ex.gradVmoyBB_tdt->Transpose()); }

	 // variation de la déformation / au ddl  
     int d_epsBB_tdtTaille = d_epsBB_tdt.Taille();		      
     for (int i=1; i<= d_epsBB_tdtTaille; i++)
         *(d_epsBB_tdt(i)) = 0.5 * (*((*(ex.d_gijBB_tdt))(i)));
//     // calcul éventuelle de la dérivée seconde de la déformation
//     if (cal_derSeconde || premier_calcul)  
//	  {// recup des variations secondes de la déformation
//	   int d2_epsBB_tdtTaille1 = d2_epsBB_tdt.Taille1();
//       for (int i=1; i<= d2_epsBB_tdtTaille1; i++)
//        for (int j=1; j<= i; j++)  // symétrie du  tableau et du tenseur
//	      { *(d2_epsBB_tdt(i,j)) = 0.5 * (*((*(ex.d2_gijBB_tdt))(i,j))) ;
//	        *(d2_epsBB_tdt(j,i)) = *(d2_epsBB_tdt(i,j)) ; 
//	       }
//	   }

     // vérification éventuelle
 //    VerifCal_implicit( ex);
     return ex;
  };
  
  	  // cas explicite : tous les parametres sont de resultats	
const Met_abstraite::Expli& Deformation::Cal_explicit_Almansi (bool 
           ,TenseurBB & epsBB_t,Tableau <TenseurBB *> & d_epsBB,TenseurBB& DepsBB
           ,TenseurBB& delta_epsBB,bool premier_calcul,const Met_abstraite::Expli& ex)
  {  
     // récup des infos sur le temps 
     const VariablesTemps& vartemps = ParaGlob::Variables_de_temps();
     double deltat=vartemps.IncreTempsCourant();
     double temps =vartemps.TempsCourant();
// modif du 3 nov 2014 -> introduction de la sauvegarde de D : D_BB_t et D_BB_tdt
//     // dans le cas où le pas de temps est trop petit on met un message d'erreur
//     if (Dabs(deltat) < ConstMath::trespetit)
//     	{ cout << "\n erreur le pas de temps " << deltat
//     	       << " est trop faible pour un calcul correcte de la vitesse moyenne";
//     	  if (ParaGlob::NiveauImpression() > 3)
//     	       cout << "\n Deformation::Cal_explicit_Almansi( .... ";
//     	  Sortie(1);
//     	};
	 
     // dans le cas où ce n'est pas le premier calcul on récupère les datas qui ne sont pas recalculée
     if (!premier_calcul)
      { Deformation::Stmet& a_0 = saveDefResul->meti_00;  // pour simplifier
        metrique->Recup_grandeur_0(*a_0.giB_,*a_0.giH_,*a_0.gijBB_,*a_0.gijHH_,*a_0.jacobien_);
      };
     // epsBB_t  est dimensionne dans l'element              
     epsBB_t = 0.5 * ( *(ex.gijBB_t) - *(ex.gijBB_0));
	 
     if (Dabs(temps) > ConstMath::trespetit)
     // on peut utiliser le temps en dénominateur
       {
        // -- calcul de l'acroissement de déformation et de la vitesse de déformation
        // ici on ne connait que ce qui est à 0 et t, donc la procédure est un peu particulière
        // normalement on considère que le delta_epsBB est proportionnel à epsBB_t ! c-a-d le delta de 0 à t
        if (Dabs(deltat) >= ConstMath::trespetit)
          { delta_epsBB = (deltat/temps) * epsBB_t;
            DepsBB = delta_epsBB/deltat;}
        else// dans le cas où l'incrément de temps est nul on garde l'incrément total pour deltaeps
          { delta_epsBB =  epsBB_t;
            // pour la vitesse de déformation, on utilise une approximation linéaire sur le temps complet
            if (Dabs(temps) >= ConstMath::trespetit)
             { DepsBB = delta_epsBB/temps;}
            else // dans le cas où l'incrément de temps est nul on utilise le gradient de vitesse moyen
             // pour éviter une division par 0
   // non ça ne marche pas car le calcul du gradient moyen utilise lui aussi le delta t, mais la question qu'il faut
   // se poser c'est pourquoi delta t est si petit plutôt que d'essayer d'y palier
   // normalement passe jamais ici vu le test plus haut !!	      
             {
             
   // modif du 3 nov 2014 -> introduction de la sauvegarde de D : D_BB_t et D_BB_tdt
   //          DepsBB = 0.5 * ((*ex.gradVmoyBB_t) + ex.gradVmoyBB_t->Transpose());
               // on utilise la vitesse précédente
               DepsBB = (*saveDefResul->D_BB_t);
             }
          };
       }
     else // sinon le temps est nul
       { delta_epsBB.Inita(0.); // on met à 0 l'acroissement de déformation de t à tdt
         // on utilise la vitesse précédente pour la vitesse de déformation
         DepsBB = (*saveDefResul->D_BB_t);
       };
   
	 // variation de la déformation / au ddl  
     int  d_epsBBTaille = d_epsBB.Taille();                                          
     for (int i=1; i<= d_epsBBTaille; i++)
         *(d_epsBB(i)) = 0.5 * (*((*(ex.d_gijBB_t))(i))); 

     return ex;
  };
  
  	  // cas explicite à t et tdt: tous les parametres sont de resultats	
const Met_abstraite::Expli_t_tdt& Deformation::Cal_explicit_Almansi_tdt (bool 
           ,TenseurBB & epsBB_tdt,Tableau <TenseurBB *> & d_epsBB,TenseurBB& DepsBB
           ,TenseurBB& delta_epsBB,bool premier_calcul,const Met_abstraite::Expli_t_tdt& ex)
  {
     const VariablesTemps& vartemps = ParaGlob::Variables_de_temps();
     double deltat=vartemps.IncreTempsCourant();
// modif du 3 nov 2014 -> introduction de la sauvegarde de D : D_BB_t et D_BB_tdt
//     // dans le cas où le pas de temps est trop petit on met un message d'erreur
//     if (Dabs(deltat) < ConstMath::trespetit)
//     	{ cout << "\n erreur le pas de temps " << deltat
//     	       << " est trop faible pour un calcul correcte de la vitesse moyenne";
//     	  if (ParaGlob::NiveauImpression() > 3)
//     	       cout << "\n Deformation::Cal_explicit_Almansi_tdt( .... ";
//     	  Sortie(1);
//     	};
     // dans le cas où ce n'est pas le premier calcul on récupère les datas qui ne sont pas recalculée
     if (!premier_calcul)
      { Deformation::Stmet& a_0 = saveDefResul->meti_00;  // pour simplifier
        Deformation::Stmet& a_t = saveDefResul->meti_t;  // pour simplifier     
        metrique->Recup_grandeur_0_t(*a_0.giB_,*a_0.giH_,*a_t.giB_,*a_t.giH_
                                     ,*a_0.gijBB_,*a_0.gijHH_,*a_t.gijBB_,*a_t.gijHH_
                                     ,*a_t.gradVmoyBB_,*a_t.jacobien_,*a_0.jacobien_);
      };
	 // -- calcul de la déformation d'Almansi et de l'incrément de déformation	  	
     //     epsBB_tdt et delta_epsBB est dimensionne auparavant             
     epsBB_tdt = 0.5 * (*(ex.gijBB_tdt) - *(ex.gijBB_0));
     delta_epsBB = 0.5 * (*(ex.gijBB_tdt) - *(ex.gijBB_t));
     // -- calcul de la vitesse de déformation: a priori on utilise 
     // la variation du tenseur déformation sur delta t ce qui correspond à l'utilisation du gradient
     // de vitesse de déformation moyen calculé à t+(delta t)/2
     if (Dabs(deltat) >= ConstMath::trespetit)
        { DepsBB = delta_epsBB/deltat; }
// modif du 3 nov 2014 -> introduction de la sauvegarde de D : D_BB_t et D_BB_tdt
//     else // dans le cas où l'incrément de temps est nul on utilise le gradient de vitesse moyen
//          // pour éviter une division par 0
//          // sauf que pour l'instant il n'y a pas de calcul du gradient moyen dans la métrique
//          // mais il y a eu un test au début de la méthode qui fait que normalement on n'arrive pas ici !
//       { DepsBB = 0.5 * ((*ex.gradVmoyBB_tdt) + ex.gradVmoyBB_tdt->Transpose()); }
     else
        // on utilise la vitesse précédente pour la vitesse de déformation
        { DepsBB = (*saveDefResul->D_BB_t);};

     // variation de la déformation / au ddl
     int d_epsBB_tdtTaille = d_epsBB.Taille();		      
     for (int i=1; i<= d_epsBB_tdtTaille; i++)
         *(d_epsBB(i)) = 0.5 * (*((*(ex.d_gijBB_tdt))(i)));
     return ex;
  };

    // calcul de la deformation au temps donné  dans le cas d'une déformation d'Almansi
    // on suppose que les métriques en HH ou BB sont définis à 0 et au temps
void Deformation::Cal_Almansi_auTemps (Enum_dure temps, TenseurBB & epsBB)
  {  switch (temps)
      { case TEMPS_t:
         { // récup des infos
           const Met_abstraite::InfoExp_t& ex = metrique->Recup_InfoExp_t();
           // -- calcul de la déformation d'Almansi
           epsBB = 0.5 * ((*(ex.gijBB_t)) - (*(ex.gijBB_0)));
           break;
           }
        case TEMPS_tdt:
         { // récup des infos
           const Met_abstraite::InfoExp_tdt& ex = metrique->Recup_InfoExp_tdt();
           // -- calcul de la déformation d'Almansi
           epsBB = 0.5 * ((*(ex.gijBB_tdt)) - (*(ex.gijBB_0)));
           break;
           }
        case TEMPS_0:
         { // pb on ne peut pas calculer
           cout << "\n erreur: on demande la def a t=0 , elle n'est pas calculee ! ";
           if (ParaGlob::NiveauImpression() > 0)
             cout << "\n Deformation::Cal_Almansi_auTemps (.. ";
           cout << endl;
           Sortie(1);
           break;
           }
       };
  };