C CCGAHU    SOURCE    GOUNAND   07/07/05    21:15:08     5784           
      SUBROUTINE CCGAHU(LCOF,NOMLOI,
     $     FC,
     $     IMPR,IRET)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C***********************************************************************
C NOM         : CCGAHU
C DESCRIPTION : Calcul de la loi de comportement aux points de Gauss :
C               DEDU ADAP avec metrique
C               On implémente ici la première partie de la fonctionnelle
C               de Huang
C     
C LANGAGE     : ESOPE
C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
C               mél : gounand@semt2.smts.cea.fr
C***********************************************************************
C APPELES          :
C APPELE PAR       :
C***********************************************************************
C ENTREES            :
C ENTREES/SORTIES    :
C SORTIES            : -
C TRAVAIL            :
C***********************************************************************
C VERSION    : v1, 30/04/07, version initiale
C HISTORIQUE : v1, 30/04/07, création
C HISTORIQUE :
C***********************************************************************
C Prière de PRENDRE LE TEMPS de compléter les commentaires
C en cas de modification de ce sous-programme afin de faciliter
C la maintenance !
C***********************************************************************
-INC CCOPTIO
-INC CCREEL
CBEGININCLUDE SMCHAEL
      SEGMENT MCHAEL
      POINTEUR IMACHE(N1).MELEME
      POINTEUR ICHEVA(N1).MCHEVA
      ENDSEGMENT
      SEGMENT MCHEVA
      REAL*8 VELCHE(NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM)
      ENDSEGMENT
      SEGMENT LCHEVA
      POINTEUR LISCHE(NBCHE).MCHEVA
      ENDSEGMENT
CENDINCLUDE SMCHAEL
      INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM,N1
      POINTEUR FC.MCHEVA
      POINTEUR LCOF.LCHEVA
      POINTEUR JMAJAC.MCHEVA
      POINTEUR JMIJAC.MCHEVA
      POINTEUR JDTJAC.MCHEVA
      POINTEUR JMAREG.MCHEVA
      POINTEUR JMET.MCHEVA
      POINTEUR JTHE.MCHEVA
      POINTEUR JGAM.MCHEVA
      CHARACTER*8 NOMLOI
      INTEGER ICOF
*     
-INC TMXMAT
* Objets temporaires
      POINTEUR JAC.MXMAT,JT.MXMAT
      POINTEUR G.MXMAT,IG.MXMAT,H.MXMAT,HIG.MXMAT,HM1.MXMAT
      POINTEUR MM1.MXMAT,JTM.MXMAT,MJ.MXMAT
* Objets à préconditionner (éventuellement)
      POINTEUR A.MXMAT,B.MXMAT,C.MXMAT,D.MXMAT,F.MXMAT,M.MXMAT
*     
      SEGMENT MCOF
      POINTEUR COEF(IDIM,IDIM).MCHEVA
      ENDSEGMENT
      POINTEUR MET.MCOF
*     
      LOGICAL LBID,LFIRST
      INTEGER LAXSP
      REAL*8 DEUPI,XR
      REAL*8 XL,XM,XH
*     
      INTEGER IMPR,IRET
*     
* Executable statements
*     
      IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ccgahu'
C      IF (.NOT.(IDIM.EQ.1)) THEN
C         WRITE(IOIMP,*) 'IDIM=',IDIM,' ?'
C         GOTO 9999
C      ENDIF
      NLFC=FC.VELCHE(/6)
      NPFC=FC.VELCHE(/5)
      ICOF=0 
*     
*     Récupération des coefficients de la metrique
*     
      SEGINI MET
      DO IIDIM=1,IDIM
         ICOF=ICOF+1
         JMET=LCOF.LISCHE(ICOF)
         IF (ICOF.EQ.1) THEN 
            NLJM=JMET.VELCHE(/6)
            NPJM=JMET.VELCHE(/5)
         ELSE
            NLJM2=JMET.VELCHE(/6)
            NPJM2=JMET.VELCHE(/5)
            IF (NLJM2.NE.NLJM.OR.NPJM2.NE.NPJM) THEN 
               WRITE(IOIMP,*) 'Erreur grave dims JMET'
               GOTO 9999
            ENDIF
         ENDIF
         MET.COEF(IIDIM,IIDIM)=JMET
      ENDDO
      DO IIDIM=1,IDIM
         NJ=IDIM-IIDIM
         IF (NJ.GE.1) THEN
            DO JIDIM=IIDIM+1,IDIM
               ICOF=ICOF+1
               JMET=LCOF.LISCHE(ICOF)
               NLJM2=JMET.VELCHE(/6)
               NPJM2=JMET.VELCHE(/5)
               IF (NLJM2.NE.NLJM.OR.NPJM2.NE.NPJM) THEN 
                  WRITE(IOIMP,*) 'Erreur grave dims JMET2'
                  GOTO 9999
               ENDIF
               MET.COEF(IIDIM,JIDIM)=JMET
            ENDDO
         ENDIF
      ENDDO
*     
      ICOF=ICOF+1
      JTHE=LCOF.LISCHE(ICOF)
      NLJT=JTHE.VELCHE(/6)
      NPJT=JTHE.VELCHE(/5)
      I1  =JTHE.VELCHE(/4)
      I2  =JTHE.VELCHE(/3)
      IF ((I1.NE.1).OR.(I2.NE.1))
     $     THEN
         WRITE(IOIMP,*) 'Erreur dims JTHE'
         GOTO 9999
      ENDIF
*     
      ICOF=ICOF+1
      JGAM=LCOF.LISCHE(ICOF)
      NLJG=JTHE.VELCHE(/6)
      NPJG=JTHE.VELCHE(/5)
      I1  =JTHE.VELCHE(/4)
      I2  =JTHE.VELCHE(/3)
      IF ((I1.NE.1).OR.(I2.NE.1))
     $     THEN
         WRITE(IOIMP,*) 'Erreur dims JGAM'
         GOTO 9999
      ENDIF
*     
      ICOF=ICOF+1
      JMAJAC=LCOF.LISCHE(ICOF)
      NLJA=JMAJAC.VELCHE(/6)
      NPJA=JMAJAC.VELCHE(/5)
      IREF=JMAJAC.VELCHE(/4)
      IREL=JMAJAC.VELCHE(/3)
*     
      IF (IREL.NE.IDIM) THEN
         WRITE(IOIMP,*) 'Erreur dims JMAJAC'
         GOTO 9999
      ENDIF
*     
      ICOF=ICOF+1
      ICOF=ICOF+1
      ICOF=ICOF+1
      JMAREG=LCOF.LISCHE(ICOF)
      NLJR=JMAREG.VELCHE(/6)
      NPJR=JMAREG.VELCHE(/5)
      I1  =JMAREG.VELCHE(/4)
      I2  =JMAREG.VELCHE(/3)
      IF ((NLJR.NE.1).OR.(NPJR.NE.1).OR.(I1.NE.IREF).OR.(I2.NE.IREF))
     $     THEN
         WRITE(IOIMP,*) 'Erreur dims JMAREG'
         GOTO 9999
      ENDIF
*     
* Objets temporaires et à préconditionner
*     
      LDIM1=IREL
      LDIM2=IREF
      SEGINI,JAC
      SEGINI,MJ
      SEGINI,A
      SEGINI,B
      LDIM1=IREF
      LDIM2=IREL
      SEGINI,JT
      SEGINI,JTM
      LDIM1=IREF
      LDIM2=IREF
      SEGINI,G
      SEGINI,IG
      SEGINI,H
      SEGINI,HM1
      SEGINI,HIG
      SEGINI,C
      LDIM1=IREL
      LDIM2=IREL
      SEGINI,M
      SEGINI,MM1
      SEGINI,D
      SEGINI,F
*     
* Calcul de la métrique des éléments réguliers
*     
      CALL MAMA(JMAREG.VELCHE,IREF,IREF,
     $     'JTJ     ',H.XMAT,IREF,IREF,
     $     IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
*      SEGPRT,H
      LBID=.FALSE.
      CALL GEOLI2(IREF,1,1,H.XMAT,HM1.XMAT,DETH,LBID,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
*!!!      XH=1.0D0
      XH=1.D0/(SQRT(DETH))
*!!      XH=1.D0/(DETH)
*     
      LFIRST = .FALSE.
      DO ILFC=1,NLFC
         IF (NLJM.EQ.1) THEN
            ILJM=1
         ELSE
            ILJM=ILFC
         ENDIF
         IF (NLJA.EQ.1) THEN
            ILJA=1
         ELSE
            ILJA=ILFC
         ENDIF
         IF (NLJT.EQ.1) THEN
            ILJT=1
         ELSE
            ILJT=ILFC
         ENDIF
         IF (NLJG.EQ.1) THEN
            ILJG=1
         ELSE
            ILJG=ILFC
         ENDIF
         DO IPFC=1,NPFC
            IF (NPJM.EQ.1) THEN
               IPJM=1
            ELSE
               IPJM=IPFC
            ENDIF
            IF (NPJA.EQ.1) THEN
               IPJA=1
            ELSE
               IPJA=IPFC
            ENDIF
            IF (NPJT.EQ.1) THEN
               IPJT=1
            ELSE
               IPJT=IPFC
            ENDIF
            IF (NPJG.EQ.1) THEN
               IPJG=1
            ELSE
               IPJG=IPFC
            ENDIF
*
* Récupération de l'exposant multiplié par 2
* On suppose que c'est un entier
*
            XTHE=JTHE.VELCHE(1,1,1,1,IPJT,ILJT)
            XGAM=JGAM.VELCHE(1,1,1,1,IPJG,ILJG)
            XCO1=XTHE*0.5D0
            IXP1=NINT(IREF*XGAM)
            XCO2=(1.D0-XTHE)*((SQRT(DBLE(IREF)))**IXP1)
*!!!
*!!!         XCO2=1.D0 
            IXP2=NINT(1-XGAM)
            IF (LFIRST) THEN
               WRITE(IOIMP,*) 'XTHE=',XTHE
               WRITE(IOIMP,*) 'XGAM=',XGAM
               WRITE(IOIMP,*) 'XCO1=',XCO1
               WRITE(IOIMP,*) 'XCO2=',XCO2
               WRITE(IOIMP,*) 'IXP1=',IXP1,'/2'
               WRITE(IOIMP,*) 'IXP2=',IXP2
               LFIRST=.FALSE.
            ENDIF
*     
* Copie des coefficients de la métrique
*     
            DO IIDIM=1,IDIM
               JMET=MET.COEF(IIDIM,IIDIM)
               M.XMAT(IIDIM,IIDIM)=JMET.VELCHE(1,1,1,1,IPJM,ILJM)
            ENDDO
            DO IIDIM=1,IIDIM
               NJ=IDIM-IIDIM
               IF (NJ.GE.1) THEN
                  DO JIDIM=IIDIM+1,IDIM
                     JMET=MET.COEF(IIDIM,JIDIM)
                     M.XMAT(IIDIM,JIDIM)=JMET.VELCHE(1,1,1,1,IPJM,ILJM)
                     M.XMAT(JIDIM,IIDIM)=JMET.VELCHE(1,1,1,1,IPJM,ILJM)
                  ENDDO
               ENDIF
            ENDDO
*            SEGPRT,M
* Calcul de l'inverse
*     
           LBID=.FALSE.
           CALL GEOLI2(IREL,1,1,M.XMAT,MM1.XMAT,DETM,LBID,IMPR,IRET)
           IF (IRET.NE.0) GOTO 9999
           XDM=1.D0/(SQRT(DETM))
*1           XDM=1.D0
******           XDM=SQRT(DETM)
*     
* Copie du jacobien
*     
            CALL MAMA(JMAJAC.VELCHE(1,1,1,1,IPJA,ILJA),IREL,IREF,
     $           'COPIE   ',
     $           JAC.XMAT,IREL,IREF,
     $           IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
*            SEGPRT,JAC
*     
* Calcul de la métrique G
*     
*     Calcul de Jt
            CALL MAMA(JAC.XMAT,IREL,IREF,
     $           'TRANSPOS',JT.XMAT,IREF,IREL,
     $           IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
*     Calcul de MJ
            CALL MAMAMA(M.XMAT,IREL,IREL,JAC.XMAT,IREL,IREF,
     $           'FOIS    ',MJ.XMAT,IREL,IREF,
     $           IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
*     Calcul de JTM
            CALL MAMA(MJ.XMAT,IREL,IREF,
     $           'TRANSPOS',JTM.XMAT,IREF,IREL,
     $           IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
*     Calcul de G=JtMJ
            CALL MAMAMA(JT.XMAT,IREF,IREL,MJ.XMAT,IREL,IREF,
     $           'FOIS    ',G.XMAT,IREF,IREF,
     $           IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
* Calcul de l'inverse, du déterminant et trace de l'inverse de g
            LBID=.FALSE.
            CALL GEOLI2(IREF,1,1,G.XMAT,IG.XMAT,DETG,LBID,IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
            XM=SQRT(DETG)
*!            WRITE(IOIMP,*) 'XM=',XM
* Calcul de hg-1 et de sa trace
            CALL MAMAMA(H.XMAT,IREF,IREF,IG.XMAT,IREF,IREF,
     $           'FOIS    ',HIG.XMAT,IREF,IREF,IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
            CALL MARE(HIG.XMAT,IREF,IREF,'TRACE   ',
     $           XL,IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
* Calcul de A = (J+)t = MJg-1 (pseudo-inverse)
            CALL MAMAMA(MJ.XMAT,IREL,IREF,IG.XMAT,IREF,IREF,
     $           'FOIS    ',A.XMAT,IREL,IREF,IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
* Calcul de B = (J+)t hg-1 = MJg-1hg-1
            CALL MAMAMA(A.XMAT,IREL,IREF,HIG.XMAT,IREF,IREF,
     $           'FOIS    ',B.XMAT,IREL,IREF,IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
* Calcul de C = g-1hg-1
            CALL MAMAMA(IG.XMAT,IREF,IREF,HIG.XMAT,IREF,IREF,
     $           'FOIS    ',C.XMAT,IREF,IREF,IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
* Calcul de D = MJ g-1hg-1 JtM = B JtM
            CALL MAMAMA(B.XMAT,IREL,IREF,JTM.XMAT,IREF,IREL,
     $           'FOIS    ',D.XMAT,IREL,IREL,IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
* Calcul de F = MJ g-1 JtM = A JtM
            CALL MAMAMA(A.XMAT,IREL,IREF,JTM.XMAT,IREF,IREL,
     $           'FOIS    ',F.XMAT,IREL,IREL,IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
*     
* Calcul de la fonctionnelle
*     
            IF (NOMLOI.EQ.'AHUF    ') THEN
               XL1=SQRT(XL)**IXP1
               FONC1=XL1*XM
*
               XM1=XM**IXP2
               FONC2=XM1
*
               CONTRI=(XCO1*FONC1)+(XCO2*FONC2*(XH**(IXP2-1)))
*!!!!
               CONTRI=CONTRI*XDM
*     
* Calcul du résidu
*     
            ELSEIF (NOMLOI(1:4).EQ.'AHUR') THEN
               CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET)
               IF (IRET.NE.0) GOTO 9999
               CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET)
               IF (IRET.NE.0) GOTO 9999
               K=IDIM1
               L=IDIM2
               XL1=SQRT(XL)**IXP1
               XL2=IXP1*(SQRT(XL)**(IXP1-2))
               AKL1=XL1*XM*A.XMAT(K,L)
               AKL2=-1.D0*XL2*XM*B.XMAT(K,L)
               AKL=AKL1+AKL2
*
               XM2=IXP2*(XM**IXP2)
               CKL=XM2*A.XMAT(K,L)
*
               CONTRI=(XCO1*AKL)+(XCO2*CKL*(XH**(IXP2-1)))
*!!!!
               CONTRI=CONTRI*XDM
*     
* Calcul du jacobien
*     
            ELSEIF (NOMLOI(1:4).EQ.'AHUJ') THEN
               CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET)
               IF (IRET.NE.0) GOTO 9999
               CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET)
               IF (IRET.NE.0) GOTO 9999
               CALL CH2INT(NOMLOI(7:7),IDIM3,IMPR,IRET)
               IF (IRET.NE.0) GOTO 9999
               CALL CH2INT(NOMLOI(8:8),IDIM4,IMPR,IRET)
               IF (IRET.NE.0) GOTO 9999
               I=IDIM1
               J=IDIM2
               K=IDIM3
               L=IDIM4
               XL1=SQRT(XL)**IXP1
               XL2=IXP1*(SQRT(XL)**(IXP1-2))
               XL3=IXP1*(IXP1-2)*(SQRT(XL)**(IXP1-4))
               B11=XM*XL1*A.XMAT(I,J)*A.XMAT(K,L)
               B12=-1.D0*XL2*XM*B.XMAT(I,J)*A.XMAT(K,L)
               B13=XM*XL1*M.XMAT(K,I)*IG.XMAT(J,L)
               B14=-1.D0*XM*XL1*F.XMAT(K,I)*IG.XMAT(J,L)
               B15=-1.D0*XM*XL1*A.XMAT(K,J)*A.XMAT(I,L)
               BIJKL1=B11+B12+B13+B14+B15
               B20=+1.D0*XL3*XM*B.XMAT(I,J)*B.XMAT(K,L)
               B21=-1.D0*XL2*XM*A.XMAT(I,J)*B.XMAT(K,L)
               B22=-1.D0*XL2*XM*M.XMAT(K,I)*C.XMAT(J,L)
               B23=+1.D0*XL2*XM*F.XMAT(K,I)*C.XMAT(J,L)
               B24=+1.D0*XL2*XM*A.XMAT(K,J)*B.XMAT(I,L)
               B25=+1.D0*XL2*XM*D.XMAT(K,I)*IG.XMAT(J,L)
               B26=+1.D0*XL2*XM*B.XMAT(K,J)*A.XMAT(I,L)
               BIJKL2=B20+B21+B22+B23+B24+B25+B26
               BIJKL=BIJKL1+BIJKL2
*
               XM2=IXP2*(XM**IXP2)
               XM3=IXP2*IXP2*(XM**IXP2)
               D1=XM3*A.XMAT(I,J)*A.XMAT(K,L)
               D2=XM2*M.XMAT(K,I)*IG.XMAT(J,L)
               D3=-1.D0*XM2*F.XMAT(K,I)*IG.XMAT(J,L)
               D4=-1.D0*XM2*A.XMAT(K,J)*A.XMAT(I,L)
               DIJKL=D1+D2+D3+D4
*
               CONTRI=(XCO1*BIJKL)+(XCO2*DIJKL*(XH**(IXP2-1)))
*!!!!
               CONTRI=CONTRI*XDM
*!               WRITE(IOIMP,*) 'CONTRI=',CONTRI
            ELSE
               WRITE(IOIMP,*) 'Erreur grave'
               GOTO 9999
            ENDIF
            FC.VELCHE(1,1,1,1,IPFC,ILFC)=
     $           FC.VELCHE(1,1,1,1,IPFC,ILFC)+
     $           CONTRI
         ENDDO
      ENDDO
      SEGSUP,F
      SEGSUP,D
      SEGSUP,MM1
      SEGSUP,M
      SEGSUP,C
      SEGSUP,HIG
      SEGSUP,HM1
      SEGSUP,H
      SEGSUP,IG
      SEGSUP,G
      SEGSUP,JTM
      SEGSUP,JT
      SEGSUP,B
      SEGSUP,A
      SEGSUP,MJ
      SEGSUP,JAC
      SEGSUP,MET
*     
* Normal termination
*     
      IRET=0
      RETURN
*     
* Format handling
*     
*     
* Error handling
*     
 9999 CONTINUE
      IRET=1
      WRITE(IOIMP,*) 'An error was detected in subroutine ccgahu'
      RETURN
*     
* End of subroutine CCGAHU
*     
      END