// 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 "GeomPentaCom.h"
#include <math.h>
#include "GeomSeg.h"
#include "GeomQuadrangle.h"
#include "GeomTriangle.h"
#include "MathUtil.h"

// constructeur
    // le constructeur par défaut ne doit pas être utilisé
GeomPentaCom::GeomPentaCom()
  {   cout << "\n erreur le constructeur par defaut ne doit pas être utilise " ;
        cout << "\nGeomPentaCom::GeomPentaCom()  " << endl;     
        Sortie(1);
   };
    
//  la dimension est 3, on a nbi pt d'integration, nbn noeuds et 5 faces, 9 aretes
GeomPentaCom::GeomPentaCom(int nbi , int nbn, Enum_interpol interpol) :
    ElemGeomC0(3,nbi,nbn,5,9,PENTAEDRE,interpol)
   { 
     // coordonnees des points d'integration
     // ceux-ci sont obtenues par assemblage des points de surface et ceux d'épaisseur
     int nbil,nbis;
     switch (nbi)
      { case 1 : nbil = 1; nbis = 1; break; 
        case 2 : nbil = 2; nbis = 1; break; 
        case 3 : nbil = 3; nbis = 1; break;
        case 6 : nbil = 2; nbis = 3; break;
        case 8 : nbil = 2; nbis = 4; break;
        case 9 : nbil = 3; nbis = 3; break;
        case 12 : nbil = 3; nbis = 4; break;
        case 18 : nbil = 3; nbis = 6; break;
        default :
         {cout << "\n erreur le pentaedre  de nombre de point d'integration = " << nbi
               << "\n n\'est pas implante !! ";
          cout << "\nGeomPentaCom::GeomPentaCom(int nbi) " << endl;     
          Sortie(1);
         }
       };
       
      // un segment de base avec  deux noeuds
      GeomSeg  segbase(nbil,2); 
      // un triangle de base avec 3 noeuds
      // dans le cas de 3 pt d'integ de surface, on choisit des pt d'integ près des noeuds
      // et non au milieu des arrêtes d'où nbis_tri pour ce choix 
      int nbis_tri(nbis); if (nbis == 3) nbis_tri = 1000+nbis; 
      GeomTriangle triabase(nbis_tri,3); 
      // on reconstitue les coordonnées en trois dimension
      // et on calcul les poids d'intégration
      int npi=1;
      for (int niil=1;niil<=nbil;niil++)
        for (int nis=1;nis<=nbis;nis++,npi++)
          { ptInteg(npi) = Coordonnee
                    ((triabase.CoorPtInteg(nis))(1)
                    ,(triabase.CoorPtInteg(nis))(2)
                    ,(segbase.CoorPtInteg(niil))(1));
           // les poids d'integration sont obtenus par produit de ceux de surface
           // et ceux de hauteur
           WI(npi)=triabase.Wi(nis) * segbase.Wi(niil);
          };
};

// destructeur
GeomPentaCom::~GeomPentaCom() 
  { 
   };
// constructeur de copie
GeomPentaCom::GeomPentaCom(const GeomPentaCom& a) :
  ElemGeomC0(a)
   {
   };    
   
// en fonction de coordonnees locales, retourne true si le point est a l'interieur
// de l'element, false sinon
bool GeomPentaCom::Interieur(const Coordonnee& M)
  { if ((seg(1)->Interieur(M(3))) &&(face(1)->Interieur(Coordonnee(M(1),M(2)))))
      return true;
     else
      return false; 
   };
       
// en fonction de coordonnees locales, retourne le point local P, maximum intérieur à l'élément, donc sur la frontière
// dont les coordonnées sont sur la droite GM: c-a-d GP = alpha GM, avec apha maxi et P appartenant à la frontière
// de l'élément, G étant le centre de gravité, sauf si GM est nul, dans ce cas retour de M
Coordonnee GeomPentaCom::Maxi_Coor_dans_directionGM(const Coordonnee& M)
{ // on décompose la recherche en segment et facette
  double x=M(1); double y=M(2); double z=M(3);
  // tout d'abord on calcul le point maxi dans le plan 1 2 et dans z indépendemment
  Coordonnee XY(x,y); Coordonnee Z(z);
  Coordonnee Pxy= face(1)->Maxi_Coor_dans_directionGM(XY);
  Coordonnee Pz = seg(1)->Maxi_Coor_dans_directionGM(Z);
	 Coordonnee P(0.,0.,0.); // le retour
	 
	 // G le centre de gravité: il a la cote nulle et il est centre de gravité du triangle médian
	 double untiers = 1./3.;
	 Coordonnee G(untiers,untiers,0.);
	 Coordonnee Gxy(untiers,untiers); // en 2D xy

	 // le pb est de savoir si on sort sur une des faces triangulaires, ou sur une des faces rectangulaires
	 // Si le point maxi P sort par une face quadrangulaire, cela signifie que z_P est compris entre 1 et -1
	 // --- on va donc commencer par calculer le z_P
	 // on considère les points: G (z_G=0),  Pxy dont la côte est également nulle, M, et z_P ?
	 // soit un plan vertical (normale à Oxy) passant par G,M,Pxy. Il passe également par P. On note Mxy la projection
	 // de M sur le plan Oxy. Dans le plan verticale on peut faire une régle de trois: 
	 // on a z_P/z_M = GPxy/GMxy d'où z_P = z_M * (GPxy/GMxy)
	 // NB: si ||GMxy|| = 0, là c'est sur ça ne sort pas par un quadrangle
	 
	 
	 // 1)-----
	 Coordonnee GMxy(x-untiers,y-untiers); // en 2D xy
	 int fin_traitement = false; // par défaut
	 double gmxy = GMxy.Norme();
	 if (gmxy > ConstMath::trespetit)
	  { double z_P = z * ( (Pxy-Gxy).Norme()) / gmxy;
	    if (Abs(z_P) <= 1.) // là c'est ok on sort par un rectangle
		     {P(1)=Pxy(1); P(2)=Pxy(2); P(3)=z_P;
			     fin_traitement = true;// pour indiquer que c'est ok
			    };
		 };
	 // 2)----	
	 // dans le cas où fin_traitement est toujours faux, cela signifie que le cas sortie par un quadrangle n'a pas marché
	 // donc cela signifie que le point P est sur la face sup ou inf. En fait cela dépend de z_M si c'est positif alors
	 // c'est la face sup donc z_P=1 sinon c'est la face inf donc z_P=-1
	 // on considère un plan vertical qui passe par G, M et P . On a GMxy/z_M = GPxy/z_P avec z_P connu
	 // d'où GPxy = z_P * GMxy/z_M = e (en module) d'où \vec{GPxy} = e * \vec{GMxy)/||GMxy|| = z_P / z_M * \vec{GMxy) d'où P
	 if (! fin_traitement)	
	  { if (z > ConstMath::trespetit)
       { // cas d'une sortie par la facette supérieure, z_P = 1.
         Coordonnee GPxy = GMxy / z;
         P(1) = untiers + GPxy(1); P(2) = untiers + GPxy(2); P(3) = 1.;
         fin_traitement = true;// pour indiquer que c'est ok
        }
      else if (z < -ConstMath::trespetit)
        { // cas d'une sortie par la facette inférieure, z_P = -1.
          Coordonnee GPxy = (-1./z) * GMxy ;
          P(1) = untiers + GPxy(1); P(2) = untiers + GPxy(2); P(3) = -1.;
          fin_traitement = true;// pour indiquer que c'est ok
         }
      else
       // arrivée ici, cela veut dire que z est quasiment nulle, et pourtant due au 1)---
       // - soit Abs(z_P) est > 1 , ce qui n'est pas possible car on vient de trouver z quasiment nulle
       // - soit gmxy est également quasiment nulle : dans ce cas on choisit arbitrairement P(0,0,0)
		       { P.Zero();};
		 };
		 
  #ifdef MISE_AU_POINT
        if ((P(1) + P(2)) > (1. + ConstMath::petit))
     { cout << "\n petit soucis, le point maxi calcule n'est pas correcte ";
     P.Affiche(cout); 
     cout << "\n GeomPentaCom::Maxi_Coor_dans_directionGM(.." << endl;
    }; 
  #endif
/*	 
	 
// en fait il y a un pb, car la suite à l'air correcte si on cherche l'intersection de OM avec la frontière
// mais ce qu'il nous faut c'est GM avec la frontière, il faut donc quelques modifs !!
cout << "\n ***************** GeomPentaCom::Maxi_Coor_dans_directionGM( a revoir ******" << endl;	 
	 
    if ((Dabs(x) <= ConstMath::petit) && (Dabs(y) <= ConstMath::petit))
      { // cas où x et y transformées sont quasi nulles, il ne reste plus qu'à examiner le cas de z
        if (Dabs(z) <= ConstMath::petit)
         { // cas où toutes les coordonnées sont nulles, retour de M
			  // on met les coordonnées en x, y strictement à 0. pour éviter d'être légèrement à l'extérieur
           return Coordonnee (0.,0.,z);
          }
        else // cas où l'on modifie  suivant z (le point est sur le triangle sup ou inf)
         { return Coordonnee (0.,0.,Pz(1));
          }  
       }
    else // cas où x ou/et y sont non nulles
     {// on calcul de facteur de réduction suivant xy
      double facteurxy=0.;
      if (Dabs(x) <= ConstMath::petit) // donc y est non nulle
       {// cas donc où seule y a peut-être variée
        facteurxy = Pxy(2)/y;
        }
      else if (Dabs(y) <= ConstMath::petit) // donc x est non nulle
        {// cas donc où seule x a peut-être variée
         facteurxy = Pxy(1)/x; 
         }
      else // cas où x et y sont non nulles
        {// cas donc où x et  y x ont peut-être variées
         facteurxy = MiN(Pxy(1)/x,Pxy(2)/y);
         }   
      // on regarde maintenant du coté de z 
      if (Dabs(z) <= ConstMath::petit)
        { // z ne participe pas car nulle 
          return Coordonnee (Pxy(1),Pxy(2),z);
         }
      else // cas où la modification est possible suivant z
        { double facteurz= Pz(1)/z;
          // on prend le plus petit des deux facteurs
          if (facteurxy <= facteurz)
             { // c'est le facteur dans le plan qui agit
               return Coordonnee (Pxy(1),Pxy(2),facteurxy*z);
              }
          else
             { // c'est le facteur suivant z qui agit
               return Coordonnee (facteurz*x,facteurz*y,Pz(1));
              }
          }  
      } //-- fin du cas  où x ou/et y sont non nulles  
      // cas qui ne doit jamais arriver car on se situe forcément dans un des trois secteurs
      // précédemment étudiés
      cout << "\n erreur inconnue !!"
           << "\n GeomPentaCom::Maxi_Coor_dans_directionGM(..";
      Sortie(1);  */   
      return M; // pour éviter un warning
   };