C GRQUA2    SOURCE    CHAT      05/01/13    00:21:26     5004           
      SUBROUTINE GRQUA2(MCENVO,MCENEL,IGEOLI,INFOMC)
C------------------------------------------------------------------
C PROJET        :  CASTEM 2000
C     
C NOM           :  GRQUA2
C     
C DESCRIPTION   :  Calcule les coefficients à appliquer
C                  sur la fonction dont onv eut calculer le gradient et
C                  le Hessien.
C     
C LANGAGE       :  FORTRAN 77 + ESOPE 2000
C     
C AUTEUR        :  A. BECCANTINI, DRN/DMT/SEMT/LTMF
C     
C-----------------------------------------------------------------
C     
C APPELES (Outils) : ERREUR, KRIPAD
C                  
C     
C APPELES (Calcul) : GRQUA3, GRQUA4
C     
C-----------------------------------------------------------------
C     
C ENTREES       :
C
C                MCENVO : maillage centre-voisins
C
C                MCENEL : maillage de centre correspondant
C                         à MCENVO
C
C                IGEOLI : SPG de CHPOINT de C.L.
C
C SORTIES       :
C                INFOMC : 
C                         segment qui contient les infos 
C                         pour créer le MCHAML des coefficients
C
C-----------------------------------------------------------------
C     
C HISTORIQUE (Anomalies et modifications eventuelles)
C     
C HISTORIQUE     : 23.04.99, Creation
C     
C-----------------------------------------------------------------
C     
C     
      IMPLICIT INTEGER(I-N)
-INC CCOPTIO
-INC SMCOORD
-INC SMELEME
-INC SMLENTI
-INC SMLREEL
C     
C-----------------------------------------------------------------
C
      SEGMENT LISTPO
      INTEGER LISPO1(NPOILI),INDEX1(NBCEN+1)
      ENDSEGMENT
C
C**** LISTPO: list de points centres avec leur voisins
C              et position:
C
C     LISPO1(k) t.q.
C     LISPO1(INDEX1(j)) = numero global du centre j
C     LISPO1(INDEX1(j)),...,LISPO1(INDEX1(j+1))-1 = voisins du centre 
C        LISPO1(INDEX1(j))
C
C**** Des tableuax a deux indices
C
      SEGMENT MATRIX
      REAL*8 MAT(N1,N2)
      ENDSEGMENT
C
      POINTEUR MATA.MATRIX, MATU.MATRIX, MATV.MATRIX,MATCEL.MATRIX,
     &         COEF.MATRIX
C
C**** Segment INFOMC
C
      SEGMENT INFOMC
      INTEGER NTOTSU,LAST(NLISOU,2), REST(NBCEN)
      INTEGER LISPOI(NPOILI),INDEX(NBCEN+1)
      REAL*8  COEFF(NPOILI,LDIM)
      ENDSEGMENT
C
C**** INFOMC segment qui contient les informations pour créer le
C             MCHAML pour le calcul du gradient et de l'hessien
C
C     NTOTSU = nombre total de maillages elementaires
C      
C     LAST(i,2) = nombre de elts avec i-noueds
C     LAST(i,1) = si LAST(i,2) != 0, numero local du premier noeuds
C                 avec i-1 elts; les autres sont dans REST(LAST(i,1)),
C                 REST(REST(LAST(i,1))) si != 0
C     si j <= NBCEN, INDEX(j) contient la position du centre j dedans 
C     LISPOI(k), i.e.
C     LISPOI(INDEX(j)) = numero global du centre j
C     LISPOI(INDEX(j)),...,LISPOI(INDEX(j+1))-1 = voisins du centre 
C        LISPOI(INDEX(j))
C
C     N.B. Les informations ici sont redondantes, i.e.
C     si j appartient à la liste chaine qui parte de LAST(i,1),
C     INDEX(j+1)-INDEX(j) = i
C     
C     Mais ca nous permet de definir le MCHAML avec beaucoup de
C     facilité!
C
C     LDIM = 5 en 2D, 9 en 3D
C     COEFF(j,*) = (d()/dx,d()/dy,d2()/dxdx,d2()/dydy,d2()/dxdy)_j en 2D
C     COEFF(j,*) = (d()/dx,d()/dy,d()/dz,d2()/dxdx,d2()/dydy,d2()/dzdz,
C                   d2()/dxdy,d2()/dxdz,d2()/dydz)_j
C
      POINTEUR MCENVO.MELEME,MCENEL.MELEME,KRILIM.MLENTI,KRICEN.MLENTI
      POINTEUR MLEVOI.MLENTI,MLEVO2.MLENTI,MLEVTO.MLENTI,MLRDIS.MLREEL,
     &     MLRSIG.MLREEL, MLRCEL.MLREEL 
      
C     
C**** Variables de COOPTIO
C     
C      INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
C     &     ,IOTER,   IOLEC,  IOIMP,     IOCAR,  IOACQ
C     &     ,IOPER,   IOSGB,  IOGRA,     IOSAU,  IORES
C     &     ,IECHO,   IIMPI,  IOSPI
C     &     ,IDIM
CC     &     ,MCOORD
C     &     ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
C     &     ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
C     &     ,NORINC,NORVAL,NORIND,NORVAD
C     &     ,NUCROU, IPSAUV
C
      INTEGER IGEOLI, NTSOUS, NTOTP, ICEN, ILISOU,NBELEM, NBNOEU, NTCEN2
     &       , NLISOU, LDIM, NBMAX, NPOILI, NBCEN, IPOS, IVOIS
     &       , IELEM, NOCENE, NOCEN1, IPOI
     &       , NVOIMI
     &       , IPOSNE, JG, INDER, IVECT, N1, N2, NTOTP0
C
      REAL*8 CONTOL
C
C**** D'abord, on doit dimensioner le segment INFOMC
C     Boucle sur le LISOUS de MCENVO
C
      SEGACT MCENVO
      SEGACT MCENEL
C
      NTSOUS = MCENVO.LISOUS(/1)
C
C**** MCENEL est un maillage de POI1      
C
      IF((MCENEL.LISOUS(/1) .NE. 0) .OR. (MCENEL.ITYPEL .NE. 1))THEN
C
C******* 'MCENEL must contain only POI1 elements'
C
         MOTERR(1:40) = 'Subroutine grqua2.eso. 1                '
         WRITE(IOIMP,*) MOTERR(1:40)
         GOTO 9991
      ENDIF
C
C**** On calcule le nombre total du point dedans MCENVO 
C     (les points doubles sont contés deux fois!)
C
      NTOTP = 0
      NBMAX = 0
      ICEN = 0
      IF(NTSOUS .EQ. 0) NTSOUS = 1
C
      DO ILISOU = 1, NTSOUS
         IF(NTSOUS .EQ. 1)THEN
            IPT1 = MCENVO
         ELSE
            IPT1 = MCENVO.LISOUS(ILISOU)
            SEGACT IPT1
         ENDIF
         NBNOEU = IPT1.NUM(/1)
         NBMAX = MAX(NBMAX, NBNOEU)
C
C******* NBMAX = nombre max de noeuds dans un elt
C        Son valeur va dimensioner INFOMC (NLISOU)
C
         NBELEM = IPT1.NUM(/2)
         NTOTP = NTOTP + (NBNOEU * NBELEM)
         ICEN = ICEN + NBELEM
      ENDDO
C
      NTCEN2 = MCENEL.NUM(/2)
      IF(NTCEN2 .NE. ICEN)THEN
         MOTERR(1:40) = 'Subroutine grqua2.eso. 2                '
         WRITE(IOIMP,*) MOTERR(1:40)
         GOTO 9991
      ENDIF
      
C
      NBCEN = NTCEN2
      NPOILI = NTOTP
      SEGINI LISTPO
C     
C**** KRIPADs de MCENEL et de IGEOLI
C
      IF(IGEOLI .GT. 0) CALL KRIPAD(IGEOLI,KRILIM)
      CALL KRIPAD(MCENEL,KRICEN)
C
C**** Dedans KRIPAD:
C     SEGACT KRILIM
C     SEGACT KRICEN
C
C
C**** Boucle sur les LISOUS de MCENVO
C
      ICEN = 0
      IPOS = 0
C
C**** IPOS: pour placer les noeuds dedans LISTPO.LISPO1
C           pour conter le nombre de points et dimensioner INFOMC
C
      DO ILISOU = 1, NTSOUS
         IF(NTSOUS .EQ. 1)THEN
            IPT1 = MCENVO
         ELSE
            IPT1 = MCENVO.LISOUS(ILISOU)
C
C********** IPT1 déjà activés
C
         ENDIF
         NBNOEU = IPT1.NUM(/1)
         NBELEM = IPT1.NUM(/2)
         DO IELEM = 1 , NBELEM
            NOCENE = IPT1.NUM(NBNOEU,IELEM)
            ICEN = ICEN + 1
            NOCEN1 = MCENEL.NUM(1,ICEN)
            IF(NOCENE .NE. NOCEN1)THEN
               MOTERR(1:40) = 'Subroutine grqua2.eso. 3                '
               WRITE(IOIMP,*) MOTERR(1:40)
               GOTO 9991
            ENDIF
            IVOIS = 1
            IPOS = IPOS + 1
            LISTPO.INDEX1(ICEN)=IPOS
            LISTPO.LISPO1(IPOS) = NOCENE
            DO IPOI = 1, NBNOEU - 1
               NOCEN1 = IPT1.NUM(IPOI,IELEM)
               IF(KRICEN.LECT(NOCEN1) .NE. 0)THEN
                  IPOS = IPOS + 1
                  IVOIS = IVOIS + 1
                  LISTPO.LISPO1(IPOS) = NOCEN1
               ELSEIF(IGEOLI .GT. 0)THEN
                  IF(KRILIM.LECT(NOCEN1) .NE. 0)THEN
                     IPOS = IPOS + 1
                     IVOIS = IVOIS + 1
                     LISTPO.LISPO1(IPOS) = NOCEN1
                  ENDIF
               ENDIF
            ENDDO
         ENDDO
         IF(NTSOUS .GT. 1) SEGDES IPT1
      ENDDO
      LISTPO.INDEX1(NBCEN+1)=IPOS+1
      SEGDES MCENVO
      SEGDES MCENEL
      IF(IGEOLI .GT.0) SEGSUP KRILIM
C
C**** N.B: LISTPO est utilisé pour créer "les voisins des voisins";
C          donc il faut le completer avant de créer cette liste. 
C
C**** Maintenant on va s'amuser!
C
C     1) On doit faire une boucle sur les noeuds de
C        MCENEL 
C
C     2) en 2D il faut avoir au moins 5 voisins (6 noeuds au total)
C        en 3D il faut avoir au moins 9 noeuds (NVOIMI)
C
      IF(IDIM .EQ. 2)THEN
         NVOIMI = 6
      ELSE
         NVOIMI = 10
      ENDIF
C
C**** NTOTVO < NVOIMI:
C
C     Il faut élargir le stencil 
C
C     Je prends les noeuds les plus proches
C     
C     En suit, il faut assembler A et la "pseudo-inverser" 
C     Si rank(A) < NVOIMI, alors il faut repeter l'operation 
C     d'élargissement du stencil.
C
C     S'elle n'est pas inversible: sortir et obliger l'utilisateur a
C     refaire le maillage 
C
C**** Si non (si NTOTVO >=  NVOIMI)
C
C     il faut assembler A et la "pseudo-inverser" .
C
C     Si rank(A) < NVOIMI, alors il faut élargir le stencil ...
C
C**** Tous ca, on le stoke en INFOMC
C
C     NBCEN déjà defini
C
      NPOILI = IPOS
      NLISOU = NBMAX
      LDIM = NVOIMI - 1
      SEGINI INFOMC
C
C**** Object temporaires:
C
C     MLEVOI.MLENTI : pour consider les voisins une seule fois!
C                     liste chainé du noeud et ses voisins
C
C     On dois aussi dimensioner: 
C
C     MLEVO2.MLENTI: liste de points "voisins des voisins" du
C                    centre du quel on veut le gradient
C
C                    NBMAX-1 voisins qui ont NBMAX-1 voisins - le
C                    point centre meme (NBMAX -1) * (NBMAX - 2)
C
C     MLRDIS.MLREEL: distances des points MLEVO2 du point centre
C
C     MLEVTO.MLENTI: le point centre et tous les voisins utilize
C                    pour le calcul du gradient (NBMAX x (NBMAX-1))!
C
C     Matrices pour la "pseudo-inversion" de A avec la "SVD" 
C     (Singular Values Decomposition)
C
C     MATA.MATRIX   : matrice (NBMAX x (NBMAX - 1)) x (NVOIMI-1)
C
C     MATU.MATRIX   : matrice (NBMAX x (NBMAX - 1)) x (NVOIMI-1)
C
C     MATV.MATRIX   : matrice  (NVOIMI-1) x (NVOIMI-1)
C
C     MLRSIG.MLREEL : (NVOIMI-1) valeur singuliers de la matrice
C                     MATA
C
C     MLRCEL.MLREEL : segment de travail de longuer
C                     (NBMAX x (NBMAX - 1)) 
C
C     MATCEL.MATRIX : matrice (NBMAX x (NBMAX - 1)) x (NVOIMI-1)
C                     qui contient la transpose de la 
C                     pseudo-inverse de MATA
C
C     COEF.MATRIX   : matrice (NBMAX x (NBMAX - 1)) x (NVOIMI-1)
C                     qui contient la transpose de la 
C                     pseudo-inverse de MATA
C
C     N.B.: MATA = MATU . diag(MLRSIG) . transpose(MATV)
C           
C           pseudoinverse(MATA) = 
C               MATV . inv(diag(MLRSIG)) . transpose(MATU)
C
      JG=MCOORD.XCOOR(/1)/(IDIM+1)
      SEGINI MLEVOI
C
      JG = (NBMAX - 1) * (NBMAX - 2)
      SEGINI MLEVO2
      SEGINI MLRDIS
C
      JG = (NBMAX - 1) * NBMAX
      SEGINI MLEVTO
C
      N1 = (NBMAX - 1) * NBMAX
      N2 =  NVOIMI - 1
      SEGINI MATA
      SEGINI MATCEL
C
      N1 = (NBMAX - 1) * NBMAX
      N2 =  NVOIMI - 1
      SEGINI MATU
C
      N1 = (NBMAX - 1) * NBMAX
      N2 = NVOIMI - 1
      SEGINI MATV
C
      JG = NVOIMI - 1
      SEGINI MLRSIG
      SEGINI MLRCEL
C
      N1 = (NBMAX - 1) * NBMAX
      N2 = NVOIMI-1
      SEGINI COEF
C
      IPOSNE = 0
C
C     IPOSNE = position des noueuds dedans INFOMC.LISPOI
C
      DO IELEM = 1, NBCEN
C
C******* D'abord, il faut remplir MLEVTO
C                         definir NTOTP
C
         IPOS = LISTPO.INDEX1(IELEM)
         NOCENE = LISTPO.LISPO1(IPOS)
         MLEVTO.LECT(1)=NOCENE
         NTOTP = LISTPO.INDEX1(IELEM+1) - IPOS
         DO IVOIS = 1, NTOTP - 1, 1
            NOCEN1 = LISTPO.LISPO1(IPOS+IVOIS)
            MLEVTO.LECT(1+IVOIS)=NOCEN1
         ENDDO
         IF(NTOTP .LT. NVOIMI)THEN
            INDER = 1
         ELSE
            CONTOL = 1.0D-4
            CALL GRQUA3(CONTOL,NBMAX,NTOTP,NVOIMI,MLEVTO,MATA,MATU,MATV,
     &     MATCEL,MLRSIG,MLRCEL,COEF,INDER)
         ENDIF
         IF(INDER .EQ. 2)THEN
C
C************* Erreur anormale
C
            MOTERR(1:40) = 'Subroutine grqua2.eso. 4                '
            WRITE(IOIMP,*) MOTERR(1:40)
            GOTO 9991
         ELSEIF(INDER .EQ. 1)THEN
C
            NTOTP0 = NTOTP
            CONTOL = 1.0D-12
            CALL GRQUA4(LISTPO,MLEVOI,MLEVO2,MLRDIS,
     &           KRICEN,CONTOL,NBMAX,NTOTP,NVOIMI,MLEVTO,
     &           MATA,MATU,MATV,MATCEL,MLRSIG,MLRCEL,COEF,INDER)
C
            IF(INDER .EQ. 1)THEN
               GOTO 9992
            ELSEIF(INDER .EQ. 2)THEN
               GOTO 9991
            ENDIF
C     
            IF(NTOTP .GT. (NBMAX - 1) * NBMAX)THEN
C
C************* Erreur anormale
C
               MOTERR(1:40) = 'Subroutine grqua2.eso. 5                '
               WRITE(IOIMP,*) MOTERR(1:40)
               GOTO 9991
            ENDIF
C
C************* N.B.: Il faut ajuster  INFOMC
C
            NLISOU = MAX(NLISOU, NTOTP)
            NPOILI = NPOILI + NTOTP - NTOTP0
            SEGADJ INFOMC
         ENDIF
         IPOSNE = IPOSNE + 1
         INFOMC.INDEX(IELEM)=IPOSNE
         INFOMC.LISPOI(IPOSNE)=NOCENE
         DO IVECT = 1, NVOIMI - 1 , 1
            INFOMC.COEFF(IPOSNE,IVECT) = COEF.MAT(1,IVECT)
         ENDDO
         DO IVOIS= 2, NTOTP,1
            IPOSNE = IPOSNE + 1
            NOCEN1 = MLEVTO.LECT(IVOIS)
            INFOMC.LISPOI(IPOSNE)=NOCEN1
C
C********** Cette boucle est au contraire.
C           Voir s'on peut faire mieux!
C
            DO IVECT = 1, NVOIMI - 1 , 1
               INFOMC.COEFF(IPOSNE,IVECT) = COEF.MAT(IVOIS,IVECT)
            ENDDO
         ENDDO
C     
         IF(INFOMC.LAST(NTOTP,2) .EQ. 0) 
     $        INFOMC.NTOTSU = INFOMC.NTOTSU + 1
         INFOMC.LAST(NTOTP,2) = INFOMC.LAST(NTOTP,2) + 1
         INFOMC.REST(IELEM) = INFOMC.LAST(NTOTP,1)
C
C******* Si INFOMC.LAST(NTOTP,2)=0, alors INFOMC.LAST(NTOTP,1)=0
C
         INFOMC.LAST(NTOTP,1) = IELEM         
      ENDDO
      INFOMC.INDEX(NBCEN+1)=IPOSNE+1
      SEGSUP LISTPO
      SEGDES INFOMC
      SEGSUP KRICEN
      SEGSUP MLEVOI
      SEGSUP MLEVO2
      SEGSUP MLRDIS
      SEGSUP MLEVTO
      SEGSUP MATA
      SEGSUP MATCEL
      SEGSUP MATU
      SEGSUP MATV
      SEGSUP MLRSIG
      SEGSUP MLRCEL
      SEGSUP COEF
C     
 9999 RETURN
C
C**** Erreurs
C
 9991 CONTINUE
C
C**** Message d'erreur standard
C     5 3
C     Erreur anormale.contactez votre support
C
      CALL ERREUR(5)
      RETURN
C
C**** Donnés incompatibles
C
 9992 CONTINUE
      MOTERR(1:40) = 'Le maillage réalisé ne permet pas une   '
      WRITE(IOIMP,*) MOTERR(1:40)
      MOTERR(1:40) = 'reconstruction quadratique exacte.      '
      WRITE(IOIMP,*) MOTERR(1:40)
C
C**** Message d'erreur standard
C     21 2
C     Données incompatibles
C
      CALL ERREUR(21)
      RETURN
C
C**** Option non implementée
C
 9993 CONTINUE
C
C**** Message d'erreur standard
C     251 2
C     Tentative d'utilisation d'une option non implémentée
C
      CALL ERREUR(251)
      RETURN
C
      END