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

// constructeur
    // le constructeur par défaut ne doit pas être utilisé
GeomTetraCom::GeomTetraCom()
  {   cout << "\n erreur le constructeur par défaut ne doit pas être utilisé " ;
        cout << "\nGeomTetraCom::GeomTetraCom()  " << endl;     
        Sortie(1);
   };
    // il y a 1 points d'integration par défaut 
    // nbi : le nombre de point d'intégration voulu 
// dimension de l'espace : 3, 4 faces, 6 arêtes        
GeomTetraCom::GeomTetraCom(int nbi, int nbe, Enum_interpol interpol) :
    ElemGeomC0(3,nbi,nbe,4,6,TETRAEDRE,interpol)
{// coordonnees des points et poids d'integration
	switch (nbi)
	  {case 1 :
     { // 1 point : (ordre 1)
      double a = 1./4.;
      ptInteg(1) = Coordonnee(a,a,a);
      WI(1)= 1./6.;
	     break;
     }
    case 4 :
     { // 4 points : (ordre 2) 
      double a = (5. - sqrt(5.))/20.; double b = (5.+3.*sqrt(5.))/20.;
      ptInteg(1) = Coordonnee(a,a,a);
      ptInteg(2) = Coordonnee(a,a,b);
      ptInteg(3) = Coordonnee(a,b,a);
      ptInteg(4) = Coordonnee(b,a,a);
      for (int i =1;i<=nbi;i++)
        WI(i)= 1./24.;
	     break;
     }
    case 5 :
     { // 5 points : (ordre 3)  
      double a = 1./4.; double b=1./6.; double c = 1./2.;
      ptInteg(1) = Coordonnee(a,a,a);
      ptInteg(2) = Coordonnee(b,b,b);
      ptInteg(3) = Coordonnee(b,b,c);
      ptInteg(4) = Coordonnee(b,c,b);
      ptInteg(5) = Coordonnee(c,b,b);
      WI(1)= -2./15.;
      for (int i =2;i<=Nbi();i++)
        WI(i)= 3./40.;
	     break;
	    }
    case 15 : // formule page 371 de Dhatt et Touzot en tenant compte d'une modif (dans les i !),
     { // 15 points : (ordre 5) 
      double a = 1./4.; double b1=(7.+sqrt(15.))/34.;
      double b2=(7.-sqrt(15.))/34.;double c1=(13.+3.*sqrt(15.))/34.;
      double c2=(13.-3.*sqrt(15.))/34.;double d=(5.-sqrt(15.))/20.;
      double e=(5.+sqrt(15.))/20.; 
      ptInteg(1) = Coordonnee(a,a,a);
      ptInteg(2) = Coordonnee(b2,b2,b2);
      ptInteg(3) = Coordonnee(c1,b2,b2);
      ptInteg(4) = Coordonnee(b2,c1,b2);
      ptInteg(5) = Coordonnee(b2,b2,c1);
      ptInteg(6) = Coordonnee(b1,b1,b1);
      ptInteg(7) = Coordonnee(c2,b1,b1);
      ptInteg(8) = Coordonnee(b1,c2,b1);
      ptInteg(9) = Coordonnee(b1,b1,c2);
      ptInteg(10) = Coordonnee(d,d,e);
      ptInteg(11) = Coordonnee(d,e,d);
      ptInteg(12) = Coordonnee(e,d,d);
      ptInteg(13) = Coordonnee(d,e,e);
      ptInteg(14) = Coordonnee(e,d,e);
      ptInteg(15) = Coordonnee(e,e,d);
      WI(1)= 8./405.;
      WI(2)= (2665.+ 14.)/226800.;
      WI(3)= (2665.- 14.)/226800.;
      for (int i =4;i<=Nbi();i++)
        WI(i)= 5./567.;
	     break;
     }
    default :
     {	cout << "\nErreur : valeur incorrecte du nombre de points d'intégration demandé ! nbi = \n"
           << nbi << '\n';
      cout << "GeomTetraCom::GeomTetraCom(int nbi, int nbe)" << endl;
      Sortie(1);
     }
   };
};


// destructeur
GeomTetraCom::~GeomTetraCom() 
  { 
   };
// constructeur de copie
GeomTetraCom::GeomTetraCom(const GeomTetraCom& a) :
  ElemGeomC0(a)
   {
   };    

  
    
// en fonction de coordonnees locales, retourne true si le point est a l'interieur
// de l'element, false sinon
bool GeomTetraCom::Interieur(const Coordonnee& M)
  {  if (Dabs(M(3)) <= 1.) 
       // tout d'abord il faut que la troisième coordonnée soit comprise entre 0 et 1
       // si oui le triangle horizontale équivalent a des maxi qui vont de 0 à 0.- M(3)
       {if (Dabs(M(2)) <= 1.-M(3))
         // deuxième condition il faut que le y soit compris entre 0 et 1.-z
         // d'où un x qui doit être compris entre 0 et 1.-y-z
         {if (Dabs(M(1)) <= 1.-M(3)-M(2))
            return true;
         };
		     };
		   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 GeomTetraCom::Maxi_Coor_dans_directionGM(const Coordonnee& M)
  { // le principe est de faire un traitement en fonction du secteur dans lequel se trouve le point
    double x=M(1); double y=M(2); double z=M(3);
    // on appelle G le centre de gravité, A(1,0,0), B(0,1,0), C(0,0,0), D(0,0,1)
    double unquart = 0.25;
    Coordonnee G(unquart,unquart,unquart);
	 // 
    if ((x+y+2.*z-1. >=0.) && (2.*x + y +z -1 >=0.) && (x+2.*y+z-1. >= 0))
     { // cas où l'on est dans le secteur limité par les demies droites ]A,G] et ]B,G] et ]D,G]
	    // méthodologie:
	      //1) on calcul la normale au plan: GAB->proportionnelle à (1,1,2), on doit être à l'avant de ce plan: AM.N >= 0 -> x+y+2z-1 >=0 
		   //2) on calcul la normale au plan: GBD->proportionnelle à (2,1,1), on doit être à l'avant de ce plan: BM.N >= 0 -> 2x+y+z-1 >=0
			//3) on calcul la normale au plan: GDA->proportionnelle à (1,2,1), on doit être à l'avant de ce plan: DM.N >= 0 -> x+2y+z-1 >= 0
		 // si les 3 conditions sont ok => on doit prolonger la doite GM en un point M' appartenant au plan qui passe par ABD, de normale prop à (1,1,1)
		 // GM'. N  = alpha * GM . N  = GA . N = 0.25<3,-1,-1>(1;1;1) = 0.25 = alpha * (x+y+z - 3/4) d'où alpha et M'	
		 
       double unsSurAlpha =  (x+y+z-0.75)/0.25; 
		 // on utilise un sur alpha au lieu d'alpha pour éviter d'avoir une division par 0 au niveau du calcul d'alpha
       if (Dabs(unsSurAlpha) >= ConstMath::grand) 
        //si alpha est trop petit cela signifie que M est près de G,		 
        { return M;} 
       else
        { return (Coordonnee((M-G)/unsSurAlpha + G));}
      }
		
     else if ((2.*x+y+z-1. <= 0.) && (-x + z  >= 0.) && (-x + y >= 0))
      {// cas où l'on est dans le secteur limité par les demies droites ]C,G], ]B,G], ]D,G]
	    // méthodologie:
	      //1) on calcul la normale au plan: GCD->proportionnelle à (-1,1,0), on doit être à l'avant de ce plan: CM.N >= 0 -> -x + y >= 0 
		   //2) on calcul la normale au plan: GBC->proportionnelle à (-1,0,1), on doit être à l'avant de ce plan: CM.N >= 0 -> -x + z  >= 0.
			//3) on calcul la normale au plan: BGD->proportionnelle à (-2,-1,-1), on doit être à l'avant de ce plan: BM.N >= 0 -> 2x+y+z-1. <= 0.
		 // si les 3 conditions sont ok => on doit prolonger la doite GM en un point M' appartenant au plan qui passe par BCD, de normale prop à (1,0,0)
		 // GM'. N  = alpha * GM . N  = GC . N = -0.25 = alpha * (x-0.25) d'où alpha et M'	
		 
       double unsSurAlpha =  (0.25-x)/0.25; 
		 // on utilise un sur alpha au lieu d'alpha pour éviter d'avoir une division par 0 au niveau du calcul d'alpha
       if (Dabs(unsSurAlpha) >= ConstMath::grand) 
        //si alpha est trop petit cela signifie que M est près de G,		 
        { return M;} 
       else
        { return (Coordonnee((M-G)/unsSurAlpha + G));}
       }
     else if ((x-y >= 0.) && (x+2.*y+z -1. <= 0.) && (-y + z >= 0))
      {// cas où l'on est dans le secteur limité par les demies droites ]C,G], ]A,G], ]D,G]
	    // méthodologie:
	      //1) on calcul la normale au plan: CGD->proportionnelle à (1,-1,0), on doit être à l'avant de ce plan: CM.N >= 0 -> x - y >= 0 
		   //2) on calcul la normale au plan: AGD->proportionnelle à (-1,-2,-1), on doit être à l'avant de ce plan: CM.N >= 0 -> x+2.*y+z <== 0.
			//3) on calcul la normale au plan: AGC->proportionnelle à (0,-1,1), on doit être à l'avant de ce plan: BM.N >= 0 -> -y+z >= 0.
		 // si les 3 conditions sont ok => on doit prolonger la doite GM en un point M' appartenant au plan qui passe par ACD, de normale prop à (0,1,0)
		 // GM'. N  = alpha * GM . N  = GC . N = -0.25 = alpha * (y-0.25) d'où alpha et M'	
		 
       double unsSurAlpha =  (0.25-y)/0.25; 
		 // on utilise un sur alpha au lieu d'alpha pour éviter d'avoir une division par 0 au niveau du calcul d'alpha
       if (Dabs(unsSurAlpha) >= ConstMath::grand) 
        //si alpha est trop petit cela signifie que M est près de G,		 
        { return M;} 
       else
        { return (Coordonnee((M-G)/unsSurAlpha + G));}
       }

     else if ((x+y+2.*z -1. <=  0.) && (x-z >= 0.) && (y-z >= 0.))
      {// cas où l'on est dans le secteur limité par les demies droites ]A,G], ]B,G], ]C,G]
	    // méthodologie:
	      //1) on calcul la normale au plan: AGB->proportionnelle à (-1,-1,-2), on doit être à l'avant de ce plan: AM.N >= x+y+2.*z-1 <=  0.
		   //2) on calcul la normale au plan: BGC->proportionnelle à (1,0,-1), on doit être à l'avant de ce plan: CM.N >= 0 -> x-z >= 0.
			//3) on calcul la normale au plan: CGA->proportionnelle à (0,1,-1), on doit être à l'avant de ce plan: BM.N >= 0 -> y-z >= 0.
		 // si les 3 conditions sont ok => on doit prolonger la doite GM en un point M' appartenant au plan qui passe par ABC, de normale prop à (0,0,1)
		 // GM'. N  = alpha * GM . N  = GC . N = -0.25 = alpha * (z-0.25) d'où alpha et M'
		 
       double unsSurAlpha =  -(z-0.25)/0.25; 
		 // on utilise un sur alpha au lieu d'alpha pour éviter d'avoir une division par 0 au niveau du calcul d'alpha
       if (Dabs(unsSurAlpha) >= ConstMath::grand) 
        //si alpha est trop petit cela signifie que M est près de G,		 
        { return M;} 
       else
        { return (Coordonnee((M-G)/unsSurAlpha + G));}
       }

      // 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 GeomTetraCom::Maxi_Coor_dans_directionGM(..";
      Sortie(1);     
      return M; // pour éviter un warning
   };