C EPSIX     SOURCE    MB234859  25/09/08    21:15:22     12358          

C     PROCEDURE UTILISEE DANS LE CAS D'ELEMENTS XFEM
c     POUR LE CALCUL DE la deformation
C
C*********************************************************
C       PARTIE DECLARATIVE
C*********************************************************
      SUBROUTINE EPSIX (IMODEL,IREPS2,IVADEP,IVAEPS,
     &                  UZDPG,RYDPG,RXDPG,IIPDPG,IRETER)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)

-INC PPARAM
-INC CCOPTIO

-INC SMCOORD
-INC SMMODEL
-INC SMCHAML
-INC SMINTE
-INC SMELEME
-INC SMLREEL
C
      POINTEUR    MCHEX1.MCHELM,MINT1.MINTE,MINT2.MINTE

-INC TMPTVAL

      SEGMENT WRK1
       REAL*8 XE(3,NBBB)
       REAL*8 XDDL(LRE),XSTRS(LHOOK)
      ENDSEGMENT

      SEGMENT WRK2
       REAL*8 SHPWRK(6,NBNO),BGENE(LHOOK,LRE)
c       REAL*8 LV7WRK(NBENRMA2,2,6,NBNO)
       REAL*8 LV7WRK(NBENRMA2,2,6,NBBB)
      ENDSEGMENT

      SEGMENT MRACC
         INTEGER TLREEL(NBENRMA2,NBI)
      ENDSEGMENT

      SEGMENT MTRACE
        REAL*8 TRACE(NBPTEL)
      ENDSEGMENT
      DIMENSION UDPGE(3)

      PARAMETER  (NDDLMAX=30,NBNIMAX=10)
      PARAMETER  (NBENRMAX=5)
      DIMENSION  MLRE(NBENRMAX+1)

C      write (*,*) '############################'
C      write (*,*) '#####  DEBUT DE EPSIX  #####'
C      write (*,*) '############################'

C*********************************************************
c  Introduction du point autour duquel se fait le mouvement
c  de la section en defo plane generalisee
C*********************************************************
C IIPDPG > 0 (argument) si noeud/point support defini dans le modele
C NDPGE > 0 si besoin de calcul avec point support
      IF (IIPDPG.GT.0) THEN
        IF (IFOUR.EQ.-3) THEN
          NDPGE=3
          UDPGE(1)=UZDPG
          UDPGE(2)=RYDPG
          UDPGE(3)=RXDPG
          IREF=(IIPDPG-1)*(IDIM+1)
          XDPGE=XCOOR(IREF+1)
          YDPGE=XCOOR(IREF+2)
C*        ELSE IF (IFOUR.EQ.11) THEN
C*          NDPGE=2
C*          UDPGE(1)=UZDPG
C*          UDPGE(2)=RXDPG
C*          UDPGE(3)=0.D0
C*          XDPGE=0.D0
C*          YDPGE=0.D0
C*        ELSE IF (IFOUR.EQ. 7 .OR. IFOUR.EQ. 8 .OR. IFOUR.EQ. 9 .OR.
C*     &           IFOUR.EQ.10 .OR. IFOUR.EQ.14) THEN
C*          NDPGE=1
C*          UDPGE(1)=UZDPG
C*          UDPGE(2)=0.D0
C*          UDPGE(3)=0.D0
C*          XDPGE=0.D0
C*          YDPGE=0.D0
        else
          write(ioimp,*) 'EPSI4 : erreur NDPGE'
          call erreur(5)
          return
        ENDIF
      ELSE
        NDPGE=0
        UDPGE(1)=0.D0
        UDPGE(2)=0.D0
        UDPGE(3)=0.D0
        XDPGE=0.D0
        YDPGE=0.D0
      ENDIF

C*********************************************************
c
C       RECUPERATION + ACTIVATIONS + VALEURS UTILES
c
C*********************************************************

C++++ Recuperation de la geometrie ++++++++++++++++++++++
      MELEME= IMAMOD

C++++ RECUP DES INFOS EF ++++++++++++++++++++++++++++++++
c   + OBTENUES PAR ELQUOI DANS RIGI1 PENDANT PHASE 1
C     segment INFO deja actif dans RIGI1
c   + rigi1 n appelle pas elquoi, c'est modeli qui l'a fait
c     mais du coup, on na pas de segment minte
c     (car depend de si pt de g pour rigi, pour sigma....)
c     c'est + simple de rappeler elquoi ici
      MELE  = NEFMOD
      NGAU1 = INFELE(6)
      LRE   = INFELE(9)
      LHOOK = INFELE(10)
      MINT1 = INFMOD(6)
      MINT2 = INFMOD(10)
      MFR   = INFELE(13)
      IELE  = INFELE(14)
      NDDL  = INFELE(15)
      NSTRS = INFELE(16)

c   + AUTRES INFOS
C     nbre de noeuds par element
      NBNN1 = NUM(/1)
C     nbre d elements
      NBEL1 = NUM(/2)

c REM: pour se passer du dimensionnement du nbre d'enrichissement dans
c      elquoi et le realiser localement , on pourrait ecrire:
c      LRE = NDDLMAX*NBNN1
c      NDDL= NDDLMAX

C     sous decoupage et points de Gauss de l element geometrique de base
      if (MELE.eq.263.or.mele.eq.264) then
         NGAU2 = MINT2.POIGAU(/1)
      endif
C      write(*,*) 'dim de MINT2=6,',(MINT2.SHPTOT(/2)),(MINT2.SHPTOT(/3))
C      write(*,*) 'MINT2',(MINT2.QSIGAU(iou),iou=1,NGAU)

C++++ Recup des infos d enrichissement +++++++++++++++++++
c     recup du MCHEX1 d'enrichissement
      NOBMO1 = IVAMOD(/1)
      if(NOBMO1.ne.0) then
       do iobmo1=1,NOBMO1
        if((TYMODE(iobmo1)).eq.'MCHAML') then
         MCHEX1 = IVAMOD(iobmo1)
         segact,MCHEX1
          if((MCHEX1.TITCHE).eq.'ENRICHIS') then
            MCHAM1 = MCHEX1.ICHAML(1)
            segact,MCHAM1
            goto 1000
          endif
        endif
       enddo
       write(*,*) 'Le modele est vide (absence d enrichissement)'
      else
      write(*,*) 'Il n y a pas de MCHEML enrichissement dans le Modele'
      endif

 1000 continue
c     niveau d enrichissement(s) du modele (ddls std u exclus)
c     NBENR1= 0 si std, 1 si H seul,  2 si H+F1,  3 si H+F1+F2, etc...
      if(NOBMO1.ne.0) then
        NBENR1= MCHAM1.IELVAL(/1)
      else
        NBENR1 = 0
      endif
c      write(*,*) 'niveau d enrichissement(s) du modele',NBENR1

C*********************************************************
C       INITIALISATIONS...
C*********************************************************
      IRETER = 0
c
c     preparation des tables avec:

      if(NBENR1.ne.0) then
      do ienr=1,NBENR1
c        -le nombre d'inconnues de chaque sous-zone
c         determinee depuis le nombre de fonction de forme
c ienr=  1: U+H(1+1=2), 2: U+H+F1(2+4=6), 3: U+H+F1+F2(6+4=10)
         nbniJ = 2 + ((ienr-1)*4)
         MLRE(1+ienr) = IDIM*NBNN1*nbniJ
      enddo
      endif
C     Tables + longues car 1er indice correspond au fontion de forme std
      MLRE(1)    = IDIM*NBNN1*1

      if(NBENR1.lt.(NBENRMAX+1)) then
        do ienr=(NBENR1+1),(NBENRMAX)
           MLRE(1+ienr) = 0
        enddo
      endif
c
c      ...DU SEGMENT WRK1
      NBENRMA2 = NBENRMAX
      NBBB = NBNN1
      SEGINI,WRK1

c      ...DU SEGMENT WRK2
c      NBNO = NBNI
      NBNO = LRE/IDIM
      SEGINI,WRK2
C
c      ...DU SEGMENT MRACC
      NBENRMA2 = NBENRMAX
      NBI = NBNN1
      segini,MRACC
C
C     du nombre d erreurs sur le nom des composantes
      NBERR1=0

C*********************************************************
C
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  BOUCLE SUR LES ELEMENTS
C
      DO 2000 J=1,NBEL1

c      write(*,*) '==================================='
c      write(*,*) '========= EF',J,' /NBEL1 ========='

C*********************************************************
C     POUR CHAQUE ELEMENT, ...
C*********************************************************
C
C++++ ON RECUPERE LES COORDONNEES DES NOEUDS DE L ELEMENT
      CALL DOXE(XCOOR,IDIM,NBNN1,NUM,J,XE)

C++++ NBENRJ = niveau d enrichissement de l element ++++
C     =0 si EF std   =1 si U+H   =2 si U+H+F1   =3 si U+H+F1+F2
      NBENRJ=0
      if(NBENR1.ne.0) then
      do IENR=1,NBENR1
        MELVA1 = MCHAM1.IELVAL(IENR)
        segact,MELVA1
        do I=1,NBNN1
         mlree1 = MELVA1.IELCHE(I,J)
C        on en profite pour remplir MRACC table de raccourcis pour cet element
         TLREEL(IENR,I) = mlree1
         if(mlree1.ne.0)  then
            NBENRJ = max(NBENRJ,IENR)
C           et on active la listreel
            segact,mlree1
         endif
        enddo
      enddo
      endif
      if(NBENRMA2.gt.NBENR1) then
        do IENR2=(NBENR1+1),NBENRMA2
        do I=1,NBNN1
           TLREEL(IENR2,I) = 0
        enddo
        enddo
      endif

c++++ quelques indices pour dimensionner au plus juste
c     nbre total de ddl de l'élément considéré
      NLIGRD = MLRE(1+NBENRJ)
      NLIGRP = MLRE(1+NBENRJ)
c     nbre fonction de forme=((Ni_std+Ni_enrichi)*nbnoeud)=Ninconnues/idim
      NBNI = (MLRE(1+NBENRJ)) / IDIM

C      write(*,*) 'EF',J,' NBENRJ=',NBENRJ,
C    &'donc',NLIGRD,' ddls et ',NBNI,' fonctions de forme'

C++++ ON CALCULE XDDL ++++
      MPTVAL = IVADEP
      NCOSOU = IVAL(/1)
C
      INITOT = 0
C-->> BOUCLE SUR LES niveaux d'ENRICHISSEMENTS (0:U 1:A 2:BCDE1 3:BCDE2)
      DO IENR=0,NBENRJ
*nbre de fonction(s) de ce niveau d'enrichissement (=1 si std ou H, =4 pour F1,2,...)
        if(IENR .le. 1) then
         NBNIENR = 1
        else
         NBNIENR = 4
        endif
C---->> BOUCLE SUR LES fonctions de forme de ce type d'enrichissement
        DO INI=1,NBNIENR
          INITOT = INITOT + 1
C------>> BOUCLE SUR LA DIMENSION
          DO 2220 KDIM=1,IDIM
            ICOMP = (INITOT-1)*IDIM + KDIM

c         --cas ou on n'a pas trouvé assez de composantes--
            if(ICOMP.GT.NCOSOU) GOTO 2221

            MELVAL = IVAL(ICOMP)
c         --cas ou on a pas trouvé cette composante dans cette zone du
c           chpoint solution devenu mchaml  --
            if(MELVAL.eq.0)then
c             Avait on besoin de cette composante?
c             oui, si c'est une composante obligatoire
              if(IENR.eq.0) goto 2991
c             oui, si l'un des noeuds est enrichi
              do I=1,NBNN1
                if(TLREEL(IENR,I).ne.0) goto 2991
              enddo
c             non, si c'est facultatif et qu'on n'est pas enrichi -> on saute
              goto 2220
c           ->AVERTISSEMENT puis on saute
 2991         NBERR1=NBERR1+1
              if(IIMPI.lt.1) goto 2220
c               write(IOIMP,*) 'PB OPERATEUR EPSI :'
              write(IOIMP,991) ICOMP,IENR,INI,KDIM
  991         format(2X,'ABSENCE DANS LE CHPOINT DEPLACEMENT DE LA ',I3,
     $                  'ieme COMPOSANTE (enrichissement',I3,
     $                   ', fonction',I3,', direction ',I3,
     $                   ') NECESSAIRE POUR L UN DES NOEUDS SUIVANTS :')
              write(IOIMP,*)'  noeuds :',(NUM(iou,J),iou=1,NBNN1)
              goto 2220
            endif

C---------->> BOUCLE SUR LES NOEUDS
            IGMX = VELCHE(/1)
            IEMN = MIN(J,VELCHE(/2))
            IQ = ((INITOT-1)*NBNN1*IDIM) + KDIM
            DO I=1,NBNN1
c**              IQ = ((INITOT-1)*NBNN1*IDIM) + ((I-1)*IDIM) + KDIM
              XDDL(IQ) = VELCHE(MIN(I,IGMX),IEMN)
              IQ = IQ + IDIM
            ENDDO

 2220     CONTINUE
        ENDDO
      ENDDO

c   --cas normal (toutes les composantes souhaitees etaient presentes)--
      GOTO 2223

c   --cas ou on n'a pas trouvé assez de composantes--
 2221 CONTINUE
      if (IIMPI.ge.1) then
        WRITE(IOIMP,2222) J,NCOSOU,ICOMP
 2222   FORMAT(2X,'ATTENTION : ELEMENT ',I6,
     &  ' LE CHAMP DE DEPLACEMENT CONTIENT ',I3,' COMPOSANTES',
     &  ' PAR NOEUD AU LIEU DE ',I3,' ATTENDUES')
      endif
      NDDL=NCOSOU*NBNN1
      NBENRJ=IENR

 2223 CONTINUE
c   --cas ou on a une ou des erreurs--
      IF (NBERR1.gt.0.and.J.eq.NBEL1) THEN
        write(IOIMP,*) 'OPERATEUR GRAD : ABSENCE DANS LE CHPOINT ',
     &   'DEPLACEMENT DE CERTAINES INCONNUES ATTENDUES PAR LE MODELE'
      ENDIF
c
c rem: il serait probablement interessant au niveau du tems cpu
c  d'utiliser moins de pts de Gauss lorsque l element est élastique
c  On pourrait par ex. utiliser MINT2 = infell(12) qui contient
c  le segment d'integration de l'EF std (QUA4 par ex.)
*      if((NBENRJ.eq.0).and.(MINT2.ne.0)) then
*         MINTE = MINT2
*         NBPGAU= NGAU2
*      else
         MINTE = MINT1
         NBPGAU= NGAU1
*      endif
c
c       pour les def quadratiques
        ISDJC=0
        NBPTEL=NBPGAU
        IF (IREPS2.EQ.1) SEGINI MTRACE
C
C*********************************************************
C
C>>>>>>>>>>  BOUCLE SUR LES POINTS DE GAUSS
C
      DO 2500 KGAU=1,NBPGAU
C
C*********************************************************
C     Initialisation à 0
C*********************************************************

c ZERO ne serait-il pas facultatif?
      CALL ZERO(SHPWRK,6,NBNO)

      i6zz = 3
      IF (IDIM.EQ.3) i6zz = 4

      do ienr7=1,NBENRJ
       do inod7=1,NBNN1
        do i6=1,i6zz
           LV7WRK(ienr7,1,i6,inod7) = 0.D0
           LV7WRK(ienr7,2,i6,inod7) = 0.D0
        enddo
       enddo
      enddo

c*********************************************************
c     Calcul des fonction de forme std dans repere local
c*********************************************************

ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccccc
c     (et donc sur les Ni std)
      DO 2510 I=1,NBNN1

C++++ Calcul des Ni std
c     (rappel: 1:Ni  2:Ni,qsi  3:Ni,eta  avec i=1,4)
      SHPWRK(1,I) = SHPTOT(1,I,KGAU)
      SHPWRK(2,I) = SHPTOT(2,I,KGAU)
      SHPWRK(3,I) = SHPTOT(3,I,KGAU)
      IF (IDIM.EQ.3) SHPWRK(4,I) = SHPTOT(4,I,KGAU)

 2510 CONTINUE
ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc

c*********************************************************
c     Passage des fonctions de forme std dans repere global
c*********************************************************

C++++ CALCUL DES  Ni,x Ni,y (i=1,4) + CALCUL DE det(J)
      CALL JACOBI(XE,SHPWRK,IDIM,NBNN1,DJAC)
c      if(J.eq.1.and.KGAU.eq.1)
c     &write(*,*) 'Ni(i=1,4)=',(SHPWRK(1,iou),iou=1,NBNN1)

c*********************************************************
c     Si on est pas du tout enrichi on peut sauter une grosse partie
c*********************************************************
      if(NBENRJ.eq.0) goto 2999

c*********************************************************
c     Calcul des level set + leurs derivees dans repere global
c*********************************************************

ccccc BOUCLE SUR LES enrichissements ccccccccccccccccccc
      do 2520 ienr=1,NBENRJ

c       MELVA1=MCHAM1.IELVAL(IENR)
c       segact,MELVA1

ccccc  BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccc
       DO 2521 I=1,NBNN1

C++++ Le I eme noeud est-il ienr-enrichi?
        mlree1= TLREEL(IENR,I)
        if(mlree1.eq.0)  goto 2521

c       Calcul du repere local de fissure(=PSI,PHI)
c       (rappel: 1,1:psi  1,2:phi  2,1 psi,x  ...etc...)
        do 2522 inode=1,NBNN1
c         pour le H-enrichissement, on n a pas gardé PSI (inutile)
          if(ienr.ne.1) then
c          valeur de PSI au inode^ieme noeud
           xpsi1 = mlree1.PROG(inode)
c          qu on multiplie par la valeur de Ni^std au pt de G considéré
           LV7WRK(ienr,1,1,I)= LV7WRK(ienr,1,1,I)
     &      + (SHPWRK(1,inode)*xpsi1)
           LV7WRK(ienr,1,2,I)= LV7WRK(ienr,1,2,I)
     &      + (SHPWRK(2,inode)*xpsi1)
           LV7WRK(ienr,1,3,I)= LV7WRK(ienr,1,3,I)
     &      + (SHPWRK(3,inode)*xpsi1)
      IF (IDIM.EQ.3) LV7WRK(ienr,1,4,I)= LV7WRK(ienr,1,4,I)
     &      + (SHPWRK(4,inode)*xpsi1)
c          valeur de PHI au inode^ieme noeud
           xphi1 = mlree1.PROG(NBNN1+inode)
          else
           xphi1 = mlree1.PROG(inode)
          endif
          LV7WRK(ienr,2,1,I)= LV7WRK(ienr,2,1,I)
     &      + (SHPWRK(1,inode)*xphi1)
          LV7WRK(ienr,2,2,I)= LV7WRK(ienr,2,2,I)
     &      + (SHPWRK(2,inode)*xphi1)
          LV7WRK(ienr,2,3,I)= LV7WRK(ienr,2,3,I)
     &      + (SHPWRK(3,inode)*xphi1)
      IF (IDIM.EQ.3) LV7WRK(ienr,2,4,I)= LV7WRK(ienr,2,4,I)
     &      + (SHPWRK(4,inode)*xphi1)
 2522   continue

 2521  continue
ccccc  fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc


 2520 CONTINUE
ccccc fin de BOUCLE SUR LES enrichissements cccccccccccccccc

c     on a construit
C     LV7WRK(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD)

c*********************************************************
c     Ajout des fonctions de forme d enrichissement
c     + leurs derivees dans repere global
c*********************************************************
      CALL SHAPX(MELE,LV7WRK,NBENRMAX,NBENRJ,TLREEL,SHPWRK,IRET)


c     retour a la partie commune aux EF enrichis et non enrichis
 2999 continue

C*********************************************************
C       CALCUL DE B'
C*********************************************************
c ZERO ne serait-il pas facultatif?
c        call ZERO(BGENE,LHOOK,NLIGRP)
        KB=1
C       boucle sur tous les Ni
        DO 3001 II=1,NBNI

          BGENE(1,KB)   = SHPWRK(2,II)
          BGENE(2,KB+1) = SHPWRK(3,II)
          BGENE(4,KB)   = SHPWRK(3,II)
          BGENE(4,KB+1) = SHPWRK(2,II)

          IF(IDIM.EQ.3) THEN
            BGENE(3,KB+2)=SHPWRK(4,II)
            BGENE(5,KB)  =SHPWRK(4,II)
            BGENE(5,KB+2)=SHPWRK(2,II)
            BGENE(6,KB+1)=SHPWRK(4,II)
            BGENE(6,KB+2)=SHPWRK(3,II)
          ENDIF

          KB = KB + IDIM

 3001   CONTINUE
C
c      if(J.eq.5.and.KGAU.eq.1) then
c      if(KGAU.eq.1) then
c      write(*,*) 'BGENE(1,..)=',(BGENE(1,iou),iou=1,2*NBNI)
c      write(*,*) 'BGENE(2,..)=',(BGENE(2,iou),iou=1,2*NBNI)
c      write(*,*) 'BGENE(6,..)=',(BGENE(6,iou),iou=1,2*NBNI)
c      endif
c      endif

C*********************************************************
C       CALCUL DE  B.q'
C*********************************************************

C  APPEL A LA PROCEDURE DE CALCUL
        CALL BST(BGENE,XDDL,(NBNI*IDIM),LHOOK,XSTRS)
c
c       cas de la priâe en compte des termes quadratiques
        IF (IREPS2.EQ.1)
     &  CALL BST2(SHPWRK,XDDL,XE,NBNI,IFOUR,XSTRS,TRACE,KGAU,
     &            XDPGE,YDPGE,UDPGE,NIFOUR)

C*********************************************************
C       ECRITURE DES DIFFERENTES COMPOSANTES DES DEFORMATIONS
C*********************************************************
        MPTVAL = IVAEPS
        DO 7000 ICOMP=1,LHOOK
          MELVAL = IVAL(ICOMP)
          VELCHE(KGAU,J) = XSTRS(ICOMP)
 7000   CONTINUE
C     if(KGAU.eq.1) then
C      write(*,*) J,KGAU,'EPSI(..)=',(XSTRS(iou),iou=1,LHOOK)
C     endif

 2500 CONTINUE
C FIN DE BOUCLE SUR LES POINTS DE GAUSS <<<<<<<<<<<<<<

c     quelques suppressions
c     (ici element non-incompressible=> pas besoin de MTRACE (cf epsi2)
      IF (IREPS2.EQ.1) THEN
        SEGSUP MTRACE
      ENDIF

 2000 CONTINUE
C FIN DE BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<<

C*********************************************************
C       SUPPRESSION ET DESACTIVATION DE SEGMENTS
C*********************************************************
C
      SEGSUP,WRK1,WRK2
      SEGSUP,MRACC

      END