// 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 "FrontSegQuad.h"
#include "Base.h"
#include "MathUtil.h"
#include "FrontPointF.h"
#include "Met_biellette.h"
#include "MetAxisymetrique2D.h"
#include "Util.h"

//----------------------------------------------------------------
// def des donnees commune a tous les elements
//----------------------------------------------------------------
// 2 points d'integration et 3 noeuds
GeomSeg FrontSegQuad::segment(2,3);                                     
Met_abstraite * FrontSegQuad::met = NULL;
Vecteur FrontSegQuad::phi(3);
Mat_pleine FrontSegQuad::dphi(1,3);
BaseB FrontSegQuad::giB;  
BaseH FrontSegQuad::giH;  
    
    // CONSTRUCTEURS :
FrontSegQuad::FrontSegQuad () : // par defaut
  ElFrontiere()
   { cout << "\n erreur, ce constructeur ne doit pas etre utilise "
          << "\nFrontSegQuad::FrontSegQuad ()" << endl;
     Sortie(1);
   };
    // fonction du tableau des noeuds  sommets 
FrontSegQuad::FrontSegQuad ( const Tableau <Noeud *>& tab, const DdlElement& ddlElem) :
  ElFrontiere(tab,ddlElem,3),
  refP(tab(1)->Dimension()),droite (tab(1)->Dimension()),theta(1),theta_repere(1)
  ,d_T()
  { // au premier appel on construit la metrique associee
    if ( met == NULL)
      DefMetrique();
    // définition de d_T
    int nb_ddl = 3 * tab(1)->Coord0().Dimension();
    d_T.Change_taille(nb_ddl);
  };
// de copie   
FrontSegQuad::FrontSegQuad( const FrontSegQuad& a) :
  ElFrontiere(a),refP(a.refP),droite(a.droite),theta(a.theta),theta_repere(a.theta_repere)
  ,d_T(a.d_T)
  {};

    // DESTRUCTEUR :
FrontSegQuad::~FrontSegQuad ()
 {};
// surcharge de l'affectation
ElFrontiere& FrontSegQuad::operator = ( const ElFrontiere& a)
   { if (this->TypeFrontiere() == a.TypeFrontiere())
       { this->ElFrontiere::operator=(a);
         const FrontSegQuad* b = (const FrontSegQuad*) &a;
         refP = Ref(); droite = b->DR(); theta = b->Theta();
         theta_repere = b->theta_repere;
         d_T=b->d_T;
         return *this;
       }
     else
       { cout << "\n erreur d\'affectation, le deux membres non pas le meme type ";
         cout << "\n FrontSegQuad& ElFrontiere::operator = (ElFrontiere& a) " << endl;
         Sortie (1);
         return *this;
       };
   };
     
// retourne le type de l'element frontiere
string FrontSegQuad::TypeFrontiere() const 
  { return string("FrontSegQuad");};
  
// creation d'un nouvelle element frontiere du type FrontSegQuad
ElFrontiere * FrontSegQuad::NevezElemFront() const 
  { ElFrontiere * pt;
    pt = new FrontSegQuad(*this);
    return pt; 
  };
// creation d'un nouvelle element frontiere du type FrontSegQuad
// avec des donnees differentes
ElFrontiere * FrontSegQuad::NevezElemFront
  ( const Tableau <Noeud *>& tab, const DdlElement& ddlElem ) const 
  { ElFrontiere * pt;
    pt = new FrontSegQuad(tab,ddlElem);
    return pt; 
  };
// ramene et calcul les coordonnees du point de reference de l'element
Coordonnee FrontSegQuad::Ref()
  { // le point de ref en coordonnees locale
    Coordonnee A0(0.);
    // appel de la routine de metrique, 3 noeuds
    refP = met->PointM_tdt(tabNoeud,segment.Phi_point(A0));
    return refP;
  };

// ramene  une droite tangente au point de reference
// si indic = 1 -> une droite
// ces infos sont stocke et sauvegardees dans l'element 
void FrontSegQuad::TangentRef(Droite& dr, Plan& , int& indic)
  { indic = 1;
    // le point de ref en coordonnees locale
    Coordonnee A0(0.);
    BaseB giB;
    met->BaseND_tdt(tabNoeud,segment.Dphi_point(A0),segment.Phi_point(A0),giB,giH);
    droite.Change_ptref(Ref());
//// debug
//{                  double d = giB(1).Coor().Norme();
//              //    if (d <= ConstMath::trespetit)
//                  if (d <= ConstMath::petit) // *****pour test
//                   { cout << "\n***  debug FrontSegQuad::TangentRef( : la norme du vecteur est trop petite !";
//                     cout <<"\nnorme = " << d << " newtra= " << giB(1).Coor() ;
//                    }
//}
//// fin debug
    droite.Change_vect(giB(1).Coor());
    theta(1) = 0.;theta_repere=theta;
    dr = droite;
  };
    // M est un point de la  derniere droite tangente sauvegarde dans l'element 
    //  - calcul du point M1 correspondant sur la courbe , M1 est stocke
    //  _ calcul et retour de la  droite tangente  au point M1
    // si indic = 1 -> une droite
    // ces infos sont stocke et sauvegardees dans l'element 
void FrontSegQuad::Tangent(const Coordonnee& M,Coordonnee& M1, Droite& dr, Plan& , int& indic)
  { // récup des bases au point courant projeté
    met->BaseND_tdt(tabNoeud,segment.Dphi_point(theta_repere),segment.Phi_point(theta_repere),giB,giH);
    // on incremente la coordonnee curviligne
    Coordonnee M1M =  M - droite.PointDroite();
    theta(1) += M1M * giH.Coordo(1);
//// debug
//Coordonnee M1M_save(M1M);
//Coordonnee M_save(M);
//double theta1_save = theta(1);
//double theta_repere_save = theta_repere(1);
//// fin partie debug


    // dans le cas où le point est externe à l'élément, on limite le repère de calcul au point externe
    // de l'élément dans la direction de theta
    if (!(segment.Interieur(theta)))
       // calcul d'un point  extreme de l'élément dans le sens de M
       { theta_repere = segment.Maxi_Coor_dans_directionGM(theta);
		       // on recalcule le repère local
         phi = segment.Phi_point(theta_repere); // fonctions d'interpolation au point limite dans l'élément
         met->BaseND_tdt(tabNoeud,segment.Dphi_point(theta_repere),phi,giB,giH);
			      // on recalcule les coodonnées locales
         theta(1) += M1M * giH.Coordo(1);
       }
    else 
	    // sinon le point est correcte, on peut y calculer le nouveau repère
       { theta_repere = theta;
         phi = segment.Phi_point(theta_repere); // fonctions d'interpolation au point
         met->BaseND_tdt(tabNoeud,segment.Dphi_point(theta_repere),phi,giB,giH);
       };
	 
    // calcul du point correspondant au theta_i
    Coordonnee delta_theta = theta - theta_repere;
    M1 = met->PointM_tdt(tabNoeud,phi) + delta_theta(1) * giB.Coordo(1);

    droite.Change_ptref(M1);

//   //---debug
//   if ((!isfinite(theta(1))) || (isnan(theta(1))))
//      {cout << "\n***  debug FrontSegQuad::Tangent( : theta(1) infini !";
//       cout << "\n M1M_save= " << M1M_save
//            << "\n theta1_save= "<< theta1_save
//            << " theta_repere_save= "<< theta_repere_save
//            << "\n M_save= " << M_save;
//      };
//   if ((!isfinite(theta_repere(1))) || (isnan(theta_repere(1))))
//   {cout << "\n***  debug FrontSegQuad::Tangent( : theta_repere(1) infini !";
//   };
//   //---fin debug


//// debug
//{                  double d = giB(1).Coor().Norme();
//              //    if (d <= ConstMath::trespetit)
//                  if (d <= ConstMath::petit) // *****pour test
//                   { cout << "\n***  debug FrontSegQuad::Tangent( : la norme du vecteur est trop petite !";
//                     cout <<"\nnorme = " << d << " newtra= " << giB(1).Coor() ;
//                     cout << "\n tabNoeud ";
//                     for (int i=1;i<= tabNoeud.Taille();i++)
//                       cout << "\n noe i= " << i <<": " << (*tabNoeud(i));
//                     cout << "\n theta(1)= " << theta(1);
//                     cout << "\n M1M_save= " << M1M_save
//                          << "\n theta1_save= "<< theta1_save
//                          << "\n M_save= " << M_save;
//                    }
//}
//// fin debug

    droite.Change_vect(giB(1).Coor());
    // retour
    dr = droite;
    indic = 1;
  };

// ramene une autre droite tangente  genere de maniere pseudo aleatoire
// si indic = 1 -> une droite,
// ces infos sont stocke et sauvegardees dans l'element
void FrontSegQuad::AutreTangent(Droite& dr, Plan& , int& indic)
  { // on genere un nombre entre -1 et 1
    nrand++; // pour avoir un nombre different
    srand(nrand);
    theta(1) = -1. + 2.*((double) rand())/RAND_MAX;
    // dans le cas où le point est externe à l'élément, on limite le repère de calcul au point externe
    // de l'élément dans la direction de theta
    if (!(segment.Interieur(theta)))
       // calcul d'un point  extreme de l'élément dans le sens de M
       { theta_repere = segment.Maxi_Coor_dans_directionGM(theta);
		       // on recalcule le repère local
         phi = segment.Phi_point(theta_repere); // fonctions d'interpolation au point limite dans l'élément
         met->BaseND_tdt(tabNoeud,segment.Dphi_point(theta_repere),phi,giB,giH);
		     }
	   else
	    // sinon le point est correcte, on peut y calculer le nouveau repère 
       { theta_repere = theta;
         phi = segment.Phi_point(theta_repere); // fonctions d'interpolation au point
         met->BaseND_tdt(tabNoeud,segment.Dphi_point(theta_repere),phi,giB,giH);
        };

////---debug
//if ((!isfinite(theta(1))) || (isnan(theta(1))))
//   {cout << "\n***  debug FrontSegQuad::AutreTangent( : theta(1) infini !";
//   };
//if ((!isfinite(theta_repere(1))) || (isnan(theta_repere(1))))
//{cout << "\n***  debug FrontSegQuad::AutreTangent( : theta_repere(1) infini !";
//};
////---fin debug

    // calcul du point correspondant au theta_i
    Coordonnee delta_theta = theta - theta_repere;
    Coordonnee M1 = met->PointM_tdt(tabNoeud,phi) + delta_theta(1) * giB.Coordo(1);

    droite.Change_ptref(M1);
//// debug
//{                  double d = giB(1).Coor().Norme();
//              //    if (d <= ConstMath::trespetit)
//                  if (d <= ConstMath::petit) // *****pour test
//                   { cout << "\n***  debug FrontSegQuad::AutreTangentRef( : la norme du vecteur est trop petite !";
//                     cout <<"\nnorme = " << d << " newtra= " << giB(1).Coor() ;
//                    }
//}
//// fin debug

    droite.Change_vect(giB(1).Coor());
    // retour
    dr = droite;
    indic = 1;
  };
    // ramene true si le dernier point M1 est dans l'element sinon false
    // le calcul est fait à eps relatif près
bool FrontSegQuad::InSurf(const double& eps) const
  { 
    if (Dabs(theta(1)) <= 1.+eps)
        return true;
    else
      return false;  
  };
// actualise et ramene le dernier plan tangent (ou  droite tangente) calcule
// si indic = 1 -> une droite, =2 -> un plan
// ramène éventuellement la variation du vecteur normale pour un plan en 3D ou une ligne en 2D
// dans le cas d'une ligne en 3D ramène la variation du vecteur tangent: si var_normale = true, sinon ramène NULL
Tableau <Coordonnee >* FrontSegQuad::DernierTangent(Droite& dr, Plan& , int& indic,bool avec_var)
  { 
    // dans le cas où le point est externe à l'élément, on limite le repère de calcul au point externe
    // de l'élément dans la direction de theta
    if (!(segment.Interieur(theta)))
       // calcul d'un point  extreme de l'élément dans le sens de M
       { theta_repere = segment.Maxi_Coor_dans_directionGM(theta);
		       // on recalcule le repère local
         phi = segment.Phi_point(theta_repere); // fonctions d'interpolation au point limite dans l'élément
         met->BaseND_tdt(tabNoeud,segment.Dphi_point(theta_repere),phi,giB,giH);
		     }
	   else
	    // sinon le point est correcte, on peut y calculer le nouveau repère 
       { theta_repere = theta;
         phi = segment.Phi_point(theta_repere); // fonctions d'interpolation au point
         met->BaseND_tdt(tabNoeud,segment.Dphi_point(theta_repere),phi,giB,giH);
       };
//  //---debug
//  if ((!isfinite(theta(1))) || (isnan(theta(1))))
//     {cout << "\n***  debug FrontSegQuad::DernierTangent( : theta(1) infini !";
//     };
//  if ((!isfinite(theta_repere(1))) || (isnan(theta_repere(1))))
//  {cout << "\n***  debug FrontSegQuad::DernierTangent( : theta_repere(1) infini !";
//  };
//  //---fin debug

    // calcul du point correspondant au theta_i
    Coordonnee delta_theta = theta - theta_repere;
    Coordonnee M1 = met->PointM_tdt(tabNoeud,phi) + delta_theta(1) * giB.Coordo(1);

    droite.Change_ptref(M1);
//    // debug
//    {                  double d = giB(1).Coor().Norme();
//                  //    if (d <= ConstMath::trespetit)
//                      if (d <= ConstMath::petit) // *****pour test
//                       { cout << "\n***  debug FrontSegQuad::DernierTangentRef( : la norme du vecteur est trop petite !";
//                         cout <<"\nnorme = " << d << " newtra= " << giB(1).Coor() ;
//                        }
//    }
//    // fin debug
    droite.Change_vect(giB(1).Coor());
    // prépa retour
    dr = droite;
    indic = 1;
    // dans le cas où l'on veut la variation du vecteur normal
    Tableau <Coordonnee >* pt_varT = NULL; // par défaut
    if (avec_var)
     { // calcul de la variation des vecteurs de base: donc ici d'un seul vecteur
       const Tableau <BaseB>& d_giB_tdt = met->d_BaseNat_tdt(tabNoeud,segment.Dphi_point(theta_repere),phi);
       // calcul de la variation de la tangente normée
       Util::VarUnVect_coorBN( giB(1),d_giB_tdt,giB(1).Coor().Norme(),d_T);
       if (ParaGlob::Dimension() == 2)
         { // dans le cas 2D on peut calculer la normale, on le fait donc
           int taille = d_T.Taille();
           for (int i=1;i<=taille;i++)
            { Coordonnee& dT = d_T(i); // par simplicité
              double inter = dT(1);
              dT(1) = -dT(2);
              dT(2) = inter;
            };
         };
       pt_varT = &d_T;
     };
    // retour 
    return pt_varT;	 
  };
   
// calcul les fonctions d'interpolation au dernier point de projection sauvegarde
const Vecteur& FrontSegQuad::Phi()
 { return segment.Phi_point(theta);};
   
// test si la position d'un point est du bon cote ( c-a-d hors matiere) ou non
// si le point est sur la surface, ramène false
// ramene true si hors matiere, sinon false
// le test  sur a est executer uniquement dans les cas suivants :
// dimension 3D et  frontiere 2D
// dimension 3D axi et frontière 1D
// dimension 2D et frontiere 1D
// ->>> dimension 3D et frontiere 1D, pas de verif
// ->>> autre cas ne doivent pas arriver normalement !!
// retour de r = distance du point à la surface, ligne
bool FrontSegQuad::BonCote_t( const Coordonnee& a,double& r)  const  // cas ou on utilise la frontiere a t
  { if ((refP.Dimension() == 2)|| (ParaGlob::AxiSymetrie()))
     { // def des infos de la droite   tangente a t en theta
       // devrait également marcher pour l'axisymétrie en 3D
       int dim = ParaGlob::Dimension();
		     phi = segment.Phi_point(theta_repere); // fonctions d'interpolation au point
		     // calcul des bases et du plan tangent et de giH
		     met->BaseND_tdt(tabNoeud,segment.Dphi_point(theta_repere),phi,giB,giH);
       // def des infos du  plan tangent a t en theta
       Coordonnee delta_theta = theta - theta_repere;
	      Coordonnee M1 = met->PointM_t(tabNoeud,phi) + delta_theta(1) * giB.Coordo(1);
       // calcul d'un vecteur proportionnel à la normale bien oriente: c-a-d vers l'extérieur du solide
       Vecteur Nor(dim); Nor(1) =  giB(1)(2);Nor(2) = - giB(1)(1);
       r = (M1 - a).Vect() * Nor ;
       if (r < 0.)
             {return true;}
       else // le cas = 0 signifie que le point etait deja sur la facette donc pas hors matière
             {return false;};
     }
    else // cas du 3D , la normale ne veut rien dire        
      return true;
  };

bool FrontSegQuad::BonCote_tdt( const Coordonnee& a,double& r) const   // cas ou on utilise la frontiere a tdt
  {  if ((refP.Dimension() == 2)|| (ParaGlob::AxiSymetrie()))
      // on regarde si a est du bon cote de la frontiere          
      { // calcul d'un vecteur proportionnel à la normale bien oriente: c-a-d vers l'extérieur du solide
        // devrait également marcher pour l'axisymétrie en 3D
        int dim = ParaGlob::Dimension();
        Vecteur Nor(dim); Nor(1) =  droite.VecDroite()(2);Nor(2) = - droite.VecDroite()(1);
        r = (droite.PointDroite() - a).Vect() * Nor ;
        if (r < 0.)
             {return true;}
        else // le cas = 0 signifie que le point etait deja sur la facette donc pas hors matière
             {return false;};
      }
     else // cas du 3D , la normale ne veut rien dire        
       return true;
  };    
   
    
    // affichage des infos de l'elements
void FrontSegQuad::Affiche(Enum_dure temp) const 
  { cout << "\n element frontiere de type FrontSegQuad , de noeuds sommets : ";
   
    int nbn = 3;
    switch (temp)
        {case TEMPS_tdt: for  (int i =1;i<=nbn;i++)
           cout <<" noe: " << tabNoeud(i)->Num_noeud() << " "
                << tabNoeud(i)->Coord2() << ", " ; break;
         case TEMPS_t : for  (int i =1;i<=nbn;i++)
           cout <<" noe: " << tabNoeud(i)->Num_noeud() << " "
                << tabNoeud(i)->Coord1() << ", " ; break;
         case TEMPS_0 : for  (int i =1;i<=nbn;i++)
           cout <<" noe: " << tabNoeud(i)->Num_noeud() << " "
                << tabNoeud(i)->Coord0() << ", " ; break;
         default: break;
        };
   };

// creation et ramene des pointeurs sur les frontieres de l'element frontiere
// au premier appel il y a construction, ensuite on ne fait que ramener le tableau
// à moins qu'il soit effacé
Tableau <ElFrontiere*>& FrontSegQuad::Frontiere()
 {if (tabfront.Taille() == 0)
   {tabfront.Change_taille(3);
    Tableau <Noeud *> tab(1); tab(1) = tabNoeud(1);
    DdlElement ddlE(1); 
    ddlE.Change_un_ddlNoeudElement(1,ddlElem(1));
//    ddlE(1) = ddlElem(1);
    tabfront(1) = new FrontPointF(tab,ddlE);
    tab(1) = tabNoeud(2);
    ddlE.Change_un_ddlNoeudElement(1,ddlElem(2));
//    ddlE(1) = ddlElem(2);
    tabfront(2) = new FrontPointF(tab,ddlE);
    tab(1) = tabNoeud(3);
    ddlE.Change_un_ddlNoeudElement(1,ddlElem(3));
//    ddlE(1) = ddlElem(3);
    tabfront(3) = new FrontPointF(tab,ddlE);
   }
  return tabfront;
 };
  
// cas d'un élément frontière ligne:
// ramène, une longueur approximative de l'élément (toujours > 0) : calculée à l'aide
// de la ligne représentée par une suite de segments rejoignants les noeuds 
// ramène une valeur nulle, s'il n'y a pas de ligne 
double FrontSegQuad::LongueurApprox()
{ double longueur = 0.;
  if (tabNoeud(1)->ExisteCoord2())
   {longueur = (tabNoeud(1)->Coord2() - tabNoeud(2)->Coord2()).Norme()
                   + (tabNoeud(2)->Coord2() - tabNoeud(3)->Coord2()).Norme();
   }
  else if(tabNoeud(1)->ExisteCoord1())
   {longueur = (tabNoeud(1)->Coord1() - tabNoeud(2)->Coord1()).Norme()
                   + (tabNoeud(2)->Coord1() - tabNoeud(3)->Coord1()).Norme();
   } 
  else 
   {longueur = (tabNoeud(1)->Coord0() - tabNoeud(2)->Coord0()).Norme()
                   + (tabNoeud(2)->Coord0() - tabNoeud(3)->Coord0()).Norme();
   };
  return longueur;
};
  
//----- lecture écriture de restart -----
// ceci concerne uniquement les informations spécifiques 
void FrontSegQuad::Lecture_base_info_ElFrontiere_pour_projection(istream& ent)
{ string toto; 
  ent >> toto >> toto >> droite >>  toto >> theta; 	
};
void FrontSegQuad::Ecriture_base_info_ElFrontiere_pour_projection(ostream& sort)
{ sort << " FrontSQ " ;
  sort << " dr_tg " << droite <<  " theta " << theta ;
};
  
//----------- METHODES PROTEGEES : ------------------------------------
    // definition de la metrique
void FrontSegQuad::DefMetrique()
 { // dimension d'un des noeuds
   int dim_base = ((*tabNoeud(1)).Coord0()).Dimension();
   // def des variables dont on se servira
   Tableau<Enum_variable_metrique> tab(17);
   tab(1) = iM0; tab(2) = iMt; tab(3) = iMtdt ;
   tab(4)=igiB_0;tab(5)=igiB_t;tab(6)=igiB_tdt;
   tab(7)=igiH_0;tab(8)=igiH_t;tab(9)=igiH_tdt ;
   tab(10)=igijBB_0;tab(11)=igijBB_t;tab(12)=igijBB_tdt;
   tab(13)=igijHH_0;tab(14)=igijHH_t;tab(15)=igijHH_tdt ;
   tab(16) = igradVBB_tdt;tab(17) = id_giB_tdt;
   // la définition de la métrique dépend du type d'espace
   if (ParaGlob::AxiSymetrie())
   	{// cas d'un espace de travail axisymétrique
   	 met = new MetAxisymetrique2D(dim_base,ddlElem,tab,3);	
   	}
   else
    {met = new Met_biellette(dim_base,ddlElem,tab,3);};
 };