C CCGQME    SOURCE    GOUNAND   26/01/09    21:15:11     12441          
      SUBROUTINE CCGQME(LCOF,NOMLOI,
     $     FC,
     $     IMPR,IRET)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C***********************************************************************
C NOM         : CCGQME
C DESCRIPTION : Lois de comportement aux points de Gauss :
C               Qualité du maillage : alignement et isotropie
C               cf. 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, 11/05/07, version initiale
C HISTORIQUE : v1, 11/05/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 PPARAM
-INC CCOPTIO
-INC CCREEL
-INC TNLIN
*-INC 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,IH.MXMAT,H.MXMAT,HIG.MXMAT,GIH.MXMAT
      POINTEUR ME.MXMAT,MJ.MXMAT
*
      SEGMENT MCOF
      POINTEUR COEF(IDIM,IDIM).MCHEVA
      ENDSEGMENT
      POINTEUR MET.MCOF
*
      LOGICAL LBID
      INTEGER LAXSP
      REAL*8 DEUPI,XR
      REAL*8 XL,XM
      REAL*8 DETGIH(1,1),DETH(1,1)
*
      INTEGER IMPR,IRET
*
* Executable statements
*
      IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ccgqme'
C      IF (.NOT.(IDIM.EQ.1)) THEN
C         WRITE(IOIMP,*) 'IDIM=',IDIM,' ?'
C         GOTO 9999
C      ENDIF
      NLFC=FC.WELCHE(/6)
      NPFC=FC.WELCHE(/5)
      ICOF=0
      XPET=SQRT(XPETIT)
*
*     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.WELCHE(/6)
            NPJM=JMET.WELCHE(/5)
         ELSE
            NLJM2=JMET.WELCHE(/6)
            NPJM2=JMET.WELCHE(/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.WELCHE(/6)
               NPJM2=JMET.WELCHE(/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
      JMAJAC=LCOF.LISCHE(ICOF)
      NLJA=JMAJAC.WELCHE(/6)
      NPJA=JMAJAC.WELCHE(/5)
      IREF=JMAJAC.WELCHE(/4)
      IREL=JMAJAC.WELCHE(/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.WELCHE(/6)
      NPJR=JMAREG.WELCHE(/5)
      I1  =JMAREG.WELCHE(/4)
      I2  =JMAREG.WELCHE(/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
      LDIM1=IREF
      LDIM2=IREL
      SEGINI,JT
      LDIM1=IREF
      LDIM2=IREF
      SEGINI,G
      SEGINI,IH
      SEGINI,H
      SEGINI,HIG
      SEGINI,GIH
      LDIM1=IREL
      LDIM2=IREL
      SEGINI,ME
*
* Calcul de la métrique des éléments réguliers
*
      CALL MAMA(JMAREG.WELCHE,IREF,IREF,
     $     'JTJ     ',H.XMAT,IREF,IREF,
     $     IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
*      SEGPRT,H
* Calcul de l'inverse, du déterminant et trace de l'inverse de h
      CALL GEOLI2(IREF,1,1,H.XMAT,IH.XMAT,DETH,
     $     IMPR,IRET)
      IF (IRET.NE.0) THEN
         WRITE(IOIMP,*)
     $        'Metrique des elements reguliers non inversible'
         SEGPRT,JMAREG
         SEGPRT,H
         GOTO 9999
      ENDIF
*      SEGPRT,IH
*
      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
         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
*
* Copie des coefficients de la métrique
*
            DO IIDIM=1,IDIM
               JMET=MET.COEF(IIDIM,IIDIM)
               ME.XMAT(IIDIM,IIDIM)=JMET.WELCHE(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)
                     ME.XMAT(IIDIM,JIDIM)=JMET.WELCHE(1,1,1,1,IPJM,ILJM)
                     ME.XMAT(JIDIM,IIDIM)=JMET.WELCHE(1,1,1,1,IPJM,ILJM)
                  ENDDO
               ENDIF
            ENDDO
*
* Copie du jacobien
*
            CALL MAMA(JMAJAC.WELCHE(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(ME.XMAT,IREL,IREL,JAC.XMAT,IREL,IREF,
     $           'FOIS    ',MJ.XMAT,IREL,IREF,
     $           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 gh-1, de sa trace et de son déterminant, qui peut être
*     nul. La trace ne l'est pas, sauf si tous les points de l'element
*     se confondent

            CALL MAMAMA(G.XMAT,IREF,IREF,IH.XMAT,IREF,IREF,
     $           'FOIS    ',GIH.XMAT,IREF,IREF,IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
            CALL GEOLI2(IREF,1,1,GIH.XMAT,HIG.XMAT,DETGIH,
     $           IMPR,IRET)
            IF (IRET.GT.0) GOTO 9999
            XM=SQRT(MAX(DETGIH(1,1),XZERO))
            CALL MARE(GIH.XMAT,IREF,IREF,'TRACE   ',
     $           XL,IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
            IF (ABS(XL).LE.XPET) THEN
               WRITE(IOIMP,*) 'Nil jacobian matrix, check your mesh'
               GOTO 9999
            ENDIF
*
* Calcul des qualités de maillage
*
            IF (NOMLOI.EQ.'QEQU    ') THEN
               CONTRI=XM
            ELSEIF (NOMLOI.EQ.'QALI') THEN
               IF (IREF.EQ.1) THEN
                  CONTRI=1.D0
               ELSE
                  XIREF=DBLE(IREF)
                  XNUM=XL
                  XDEN=XIREF*(XM**(2.D0/XIREF))
                  XEXP=XIREF/(2.D0*(XIREF-1.D0))
                  CONTRI=(XDEN/XNUM)**XEXP
*                  CONTRI=(XNUM/XDEN)**XEXP
C                  SEGPRT,GIH
C                  WRITE(IOIMP,*) 'TRGIH =',XL
C                  WRITE(IOIMP,*) 'DETGIH=',(1.D0/DETHIG)
C                  WRITE(IOIMP,*) 'XNUM=',XNUM
C                  WRITE(IOIMP,*) 'XDEN=',XDEN
               ENDIF
            ELSE
               WRITE(IOIMP,*) 'Erreur grave'
               GOTO 9999
            ENDIF
            FC.WELCHE(1,1,1,1,IPFC,ILFC)=
     $           FC.WELCHE(1,1,1,1,IPFC,ILFC)+
     $           CONTRI
         ENDDO
      ENDDO
      SEGSUP,ME
      SEGSUP,GIH
      SEGSUP,HIG
      SEGSUP,H
      SEGSUP,IH
      SEGSUP,G
      SEGSUP,JT
      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 ccgqme'
         RETURN
*
* End of subroutine CCGQME
*
         END