C ACCRO3    SOURCE    PV090527  26/04/30    21:15:02     12529          
      SUBROUTINE ACCRO3
C========================================================================
C   Cree la matrice de liaison entre le champ u et w
C   avec une formulation forte  => On utilise les fonctions de forme
C   et on identifie le deplacement aux noeuds de la fissure
c
c   On prend en compte les enrichissements XFEM
c
C   Creation : BP, decembre 2012
C   Modifications : ...
C
C========================================================================

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

-INC PPARAM
-INC CCOPTIO
-INC CCREEL
-INC CCHAMP
-INC CCGEOME

-INC SMCOORD
-INC SMELEME
-INC SMRIGID
-INC SMLMOTS
-INC SMMODEL
-INC SMINTE
-INC SMCHAML
-INC SMLREEL

      POINTEUR    MCHEX1.MCHELM
      external shape
      DATA EPSI/1.D-9/
      DIMENSION ICOR(6),IMEL(6),imtt(10)
      DIMENSION qsi(3),XPO(3)

C  INITIALISATION DES INCONNUES obligatoires et facultatives
      PARAMETER (NOBL=3,NFAC=9)
      CHARACTER*4 DDLOBL(NOBL),DDLFAC(NFAC),MODDL,MODDL2
      CHARACTER*4 DUAOBL(NOBL),DUAFAC(NFAC)

      DATA DDLOBL/'UX  ','UY  ','UZ  '/
      DATA DDLFAC/'AX  ','AY  ','AZ  ',
     >'B1X ','B1Y ','B1Z ',
     >'B2X ','B2Y ','B2Z '/
      DATA DUAOBL/'FX  ','FY  ','FZ  '/
      DATA DUAFAC/'FAX ','FAY ','FAZ ',
     >'FB1X','FB1Y','FB1Z',
     >'FB2X','FB2Y','FB2Z'/


c-----------------------------------------------------
c     Segment de travail
c-----------------------------------------------------
      SEGMENT MTRAV
        REAL*8 XE2(3,NBNN2)
        REAL*8 SHPP1(6,NBNN1)
        REAL*8 XE1(3,NBNN1)
        INTEGER IDEJVU(NBPTS)
      ENDSEGMENT

c    tableaux comptant le nbre d EF de chaque ddl
      PARAMETER(NDDLMAX=6)
      INTEGER NELDDL(NDDLMAX)

      if(iimpi.ge.2) then
      write(ioimp,*) '-----------------------------------------------'
      write(ioimp,*) '               ENTREE dans ACCRO3'
      write(ioimp,*) '-----------------------------------------------'
      endif

c     Preliminaires
      idim1 = idim + 1
      SEGACT,MCOORD*mod

c-----------------------------------------------------
c     RECUPERATION DU MAILLAGE MASSIF
c-----------------------------------------------------
      CALL LIROBJ('MMODEL  ',IPMODL,1,IRETOU)
      CALL ACTOBJ('MMODEL  ',IPMODL,1)
      IF(IERR.NE.0) RETURN
      MMODE1=IPMODL
      segact,MMODE1
c  Récupération du nombre de zones du modèle
      N1 = MMODE1.KMODEL(/1)
      if(N1.gt.1) write(ioimp,*) 'attention 1 seule zone a ce jour!'
      IMODE1 = MMODE1.KMODEL(1)
      segdes,MMODE1
C Récupération du maillage et du numéro d'élément du modèle
      segact,IMODE1
      nele1 = IMODE1.NEFMOD
      IPT1 = IMODE1.IMAMOD
      SEGACT IPT1
C Récupération du numéro d'élément du maillage, du nombre de noeuds et d'éléments
      iele1 = IPT1.itypel
      nbnn1 = IPT1.num(/1)
      nbel1 = IPT1.num(/2)
      mele1 = IMODE1.INFELE(1)

c-----------------------------------------------------
c     RECUPERATION DU MAILLAGE INTERFACE
c-----------------------------------------------------
      IPMAI2 = 0
      CALL LIROBJ('MMODEL  ',IPMODL,0,IRETOU)
      IF(IERR.NE.0) RETURN
      if (IRETOU.EQ.1)then
        CALL ACTOBJ('MMODEL  ',IPMODL,1)
         MMODE2=IPMODL
         segact,MMODE2
         N2 = MMODE2.KMODEL(/1)
         if(N2.gt.1) write(ioimp,*) 'attention 1 seule zone a ce jour!'
         IMODE2 = MMODE2.KMODEL(1)
         segdes,MMODE2
         segact,IMODE2
         nele2 = IMODE2.NEFMOD
         IPT2 = IMODE2.IMAMOD
         SEGACT IPT2
c pour l'instant on dit que nele = iele (marche pour iele entre 2 et 26, voir bdata.eso)
         iele2 = IPT2.itypel
         nbnn2 = IPT2.num(/1)
         nbel2 = IPT2.num(/2)
         ngau2 = IMODE2.INFELE(6)
         segdes,IMODE2
      else
C Dans le cas où on a un maillage en entrée
         CALL LIROBJ('MAILLAGE',IPMAI2,1,IRETOU)
         IF(IERR.NE.0) RETURN
         IPT2 = IPMAI2
         SEGACT IPT2
c pour l'instant on dit que nele = iele (marche pour iele entre 2 et 26, voir bdata.eso)
         iele2 = IPT2.itypel
         nele2 = iele2
         if (nele2.lt.2.or.nele2.gt.26) then
         write(ioimp,*)'element geometrique different de l element fini'
            call erreur(16)
         endif
         nbnn2 = IPT2.num(/1)
         nbel2 = IPT2.num(/2)
c     SEG2
         if (nele2.EQ.2) ngau2 = 2
c     SEG3
         if (nele2.EQ.3) ngau2 = 3
c     TRI3
         if (nele2.EQ.4) ngau2 = 1
c     TRI6
         if (nele2.EQ.6) ngau2 = 4
c     QUA4
         if (nele2.EQ.8) ngau2 = 4
c     QUA8
         if (nele2.EQ.10) ngau2 = 9
      endif

      call RESHPT(ngau2,nbnn2,iele2,nele2,0,IPTR2,IRET)
      MINTE2 = IPTR2
      segact,MINTE2


c-----------------------------------------------------
c     RECHERCHE DU MCHAML ISSU MCHEX1 D ENRICHISSEMENT
c-----------------------------------------------------
      MCHAM1=0
      NBENR2=0
      segact,IMODE1
      NOBMOD = IMODE1.IVAMOD(/1)
      IF (NOBMOD.NE.0) THEN
         DO 1002 iobmo1=1,NOBMOD
            if((IMODE1.TYMODE(iobmo1)).ne.'MCHAML')   goto 1002
            MCHEX1 = IMODE1.IVAMOD(iobmo1)
            segact,MCHEX1
            if((MCHEX1.TITCHE).ne.'ENRICHIS')  goto 1003
            MCHAM1 = MCHEX1.ICHAML(1)
            segact,MCHAM1
            NBENR2 = MCHAM1.IELVAL(/1)
            do ienr2=1,NBENR2
               MELVA1=MCHAM1.IELVAL(IENR2)
               if(MELVA1.ne.0) segact,MELVA1
            enddo
 1003       continue
            segdes,MCHEX1
 1002    CONTINUE
      ENDIF

c-------------------------------------
c     INITIALISATION DES OBJETS DE TRAVAIL : MTRAV et ITYMAT
c-------------------------------------
      segini,MTRAV
*bp : NTYMAT = (U ou H ou HB1 ou HB1B2)
*     nbre de types de matrices = NRIGEL = NTYMAT * idim
c       NTYMAT = 4
      NTYMAT = 1+NBENR2

c-------------------------------------
c     INITIALISATION DU MCOORD
c-------------------------------------
      NBPTS0 = NBPTS
      NBPTS = NBPTS0 + ((nbel2 + 1)*idim*3)
      SEGADJ,MCOORD
      NBPTS = NBPTS0
*     on compte le vrai nombre de LX ajouté via NBPTS

c--------------------------------------------------------------------
c     INITIALISATION DU SEGMENT MRIGID
c--------------------------------------------------------------------
      NRIGEL = NTYMAT*idim
      segini,MRIGID
      IFORIG = IFOUR
      MTYMAT ='RIGIDITE'
c    -on prepare le meleme
      NBSOUS = 0
      NBREF = 0

      ityty=0
c on initialise la taille matrice en fonction du type de matrice
      do ity=1,NTYMAT
      do iidim=1,idim

         ityty=ityty+1
         COERIG(ityty) = 1.D0

*     dim de la matrice RE elementaire
*     = 1 (LX) + 1 (UX) +   nbnoeud (UX)
*     = 1 (LX) + 1 (AX) +   nbnoeud (AX)
*     = 1 (LX) + 1 (AX) + 2*nbnoeud (AX+B1X)
*     = 1 (LX) + 1 (AX) + 3*nbnoeud (AX+B1X+B2X)
         nbno1 = nbnn1
         nenr1 = ity
         if(nenr1.le.2) then
           NLIGRP = 1 + 1 + nbno1
         else
           NLIGRP = 1 + 1 + ((nenr1-1)*nbno1)
         endif
         NLIGRD = NLIGRP

c       -creation du MELEME
         NBNN = 2 + nbno1
         NBELEM=0
         SEGINI,MELEME
c          ITYPEL=28
         ITYPEL=22
         IRIGEL(1,ityty) = MELEME

c       -remplissage du DESCR
         SEGINI,DESCR
         IRIGEL(3,ityty) = DESCR
         iddl=0
c        remplissage des ddl LX de la fissure
         iddl=iddl+1
         LISINC(iddl)='LX'
         LISDUA(iddl)='FLX'
         NOELEP(iddl)=1
         NOELED(iddl)=1
c        remplissage des ddl UX de la fissure (ici appelé WX) []*WX=TX
         iddl=iddl+1
         if (nenr1.eq.1) then
           LISINC(iddl)=DDLOBL(iidim)
           LISDUA(iddl)=DUAOBL(iidim)
         else
           LISINC(iddl)=DDLFAC(iidim)
           LISDUA(iddl)=DUAFAC(iidim)
         endif
         NOELEP(iddl)=2
         NOELED(iddl)=2
c        remplissage des ddl de la structure
         if (nenr1.eq.1) then
           do ino1=1,nbno1
               iddl=iddl+1
               LISINC(iddl)=DDLOBL(iidim)
               LISDUA(iddl)=DUAOBL(iidim)
               NOELEP(iddl)=2+ino1
               NOELED(iddl)=2+ino1
           enddo
         else
           do ini1=1,(nenr1-1)
           do ino1=1,nbno1
               iddl=iddl+1
               LISINC(iddl)=DDLFAC(iidim+(3*(ini1-1)))
               LISDUA(iddl)=DUAFAC(iidim+(3*(ini1-1)))
               NOELEP(iddl)=2+ino1
               NOELED(iddl)=2+ino1
           enddo
           enddo
         endif
         if(iimpi.ge.3) write(ioimp,*) ityty,(LISINC(iou),iou=1,NLIGRP)
         if(iimpi.ge.3) write(ioimp,*) ityty,(NOELEP(iou),iou=1,NLIGRP)
         SEGDES,DESCR

c       -initialisation du XMATRI
         NELRIG=0
         rigrel=0
         SEGINI,XMATRI
         IRIGEL(4,ityty) = XMATRI
         IRIGEL(5,ityty) = NIFOUR
         IRIGEL(6,ityty) = 0
         IRIGEL(7,ityty) = 0
         IRIGEL(8,ityty) = 0

      enddo
      enddo


c----------------------------------------------------------------------
c     1. RECHERCHE DES ELEMENTS DE STRUCTURE CONTENANT DES POINTS DE GAUSS
c        DES ELEMENTS DE LA FISSURE
c     2. REMPLISSAGE DU MRIGID (XMATRI et MELEME)
c----------------------------------------------------------------------
      iaccro=0
      NODES=0
C
c==== Boucle sur les elements de fissure ==============================
      DO 1100 iem2=1,nbel2

         call doxe(xcoor,idim,nbnn2,ipt2.num,iem2,xe2)
         nbenrj = 0

c======= Boucle sur les noeuds de fissure ============================
         DO 1132 ino2=1,nbnn2

c   on n'attache qu'une seule fois chaque noeud
            inode2 = IPT2.NUM(ino2,iem2)
            if(IDEJVU(inode2).ne.0) goto 1132
            IDEJVU(inode2)=1
            NODES=NODES+1

c   récupération des coordonnees du point de gauss dans le repère global
            XPO(1) = xe2(1,ino2)
            XPO(2) = xe2(2,ino2)
            XPO(3) = xe2(3,ino2)

c---------- Boucle sur les elements de structure ----------------------
            DO 1131 iem1=1,nbel1

c   si pas d'enrichissement, on travaille sur tous les elements
               if(MCHAM1.eq.0) goto 1133
c   on saute les elements non enrichi car a priori ne contiennent pas la fissure
               do ienr2=1,NBENR2
                  MELVA1=MCHAM1.IELVAL(IENR2)
                  if(MELVA1.ne.0) then
                    do inode1=1,nbnn1
                  if(MELVA1.IELCHE(inode1,iem1).ne.0) goto 1133
                    enddo
                  endif
               enddo
               goto 1131
 1133          continue

c   recuperation des coordonnées des noeuds de IPT1 : xe1 (dans le repère x,y,z)
               call doxe(xcoor,idim,nbnn1,ipt1.num,iem1,xe1)

               qsi(1) = 0.D0
               qsi(2) = 0.D0
               qsi(3) = 0.D0
c   calcul des fonctions de formes de IPT1 au pt de Gauss de IPT2
               call QSIJS(xe1,iele1,nbnn1,idim,XPO,SHPP1,qsi,IRET)

c   test pour savoir si PG est dans EF de IPT1
               DO ino1=1,NBNN1
                  if (SHPP1(1,ino1).LT.-1.01D-7) then
                     go to 1131
                  endif
               ENDDO
c ON a trouvé : l'iem1 élément de structure contient ce noeud de fissure
               IDEJVU(inode2)=10
               iaccro=iaccro+1

c     DETECTION DU TYPE D'ENRICHISSEMENT MAX DE CET ELEMENT = nbenrj
               DO 3001 IENR2=1,NBENR2
                 MELVA1=MCHAM1.IELVAL(IENR2)
                 IF(MELVA1.eq.0) GOTO 3001
                 DO 3002 ino1=1,nbnn1
                   MLREEL = MELVA1.IELCHE(ino1,iem1)
c                  Test pour savoir si le noeud est enrichi
                   IF(MLREEL.eq.0) GOTO 3002
                   nbenrj=max(nbenrj,IENR2)
 3002            continue
 3001          continue

               if(iimpi.ge.3) write(ioimp,*) 'EF fissure ',iem2,
     &         ' ptdeG ',ino2,' -> EF MASSIF ',iem1,' nbenrj=',nbenrj


c Remplissage du MRIGID

c           ---Boucle sur la Partie standard et enrichie ---
               do 6000 ity=1,min(2,NTYMAT)
c           ---Boucle sur la dimension ---
c              rem : utile uniquement pour le meleme
c          (on pourrait garder le meme xmatri pour les differents iidim)
               do 6001 iidim=1,idim

                 if(ity.eq.1) then
                   ityty = iidim
                 else
                   ityty = (nbenrj*idim) + iidim
                 endif
                 MELEME = IRIGEL(1,ityty)
                 XMATRI = IRIGEL(4,ityty)
                 NBELEM = NUM(/2)+1
                 NLIGRD = RE(/1)
                 NLIGRP = RE(/2)
                 NELRIG = RE(/3)+1
                 segadj,MELEME
                 rigrel=0
                 segadj,XMATRI

c Remplissage du MELEME
c                traitement du LX
                 NBPTS = NBPTS + 1
                 if(NBPTS.gt.(XCOOR(/1)/idim1)) then
                   NBPTS0 = NBPTS
                   NBPTS = NBPTS0 + (nbel2 + 1)
                   SEGADJ,MCOORD
                   NBPTS = NBPTS0
                 endif
                 NUM(1,NBELEM) = NBPTS
                 XCOOR((NBPTS-1)*idim1 +1) = XCOOR((inode2-1)*idim1 +1)
                 XCOOR((NBPTS-1)*idim1 +2) = XCOOR((inode2-1)*idim1 +2)
                 if(idim.eq.3)
     &           XCOOR((NBPTS-1)*idim1 +3) = XCOOR((inode2-1)*idim1 +3)
c                traitement du noeud de la fissure
                 NUM(2,NBELEM) = inode2
c                traitement des noeuds de la structure
                 inono=2
cbp inutile      do j2=1,max(1,nbenrj)
cbp inutile car noelep boucle deja sur les enrichissement
                 do ino1=1,nbnn1
                   inono=inono+1
                   NUM(inono,NBELEM) = IPT1.NUM(ino1,iem1)
                 enddo
c                enddo

c Remplissage du XMATRI
c                traitement du terme LX - ddl fissure
                 RE(1,2,NELRIG)=1.d0
                 RE(2,1,NELRIG)=1.d0
c                traitement des terme LX - ddl structure
                 inono=2
c                UX = UX seulement
                 if(ity.eq.1) then
                   do ino1=1,nbnn1
                     inono=inono+1
                     RE(1,inono,NELRIG)=-1.d0*SHPP1(1,ino1)
                     RE(inono,1,NELRIG)=-1.d0*SHPP1(1,ino1)
                   enddo
                 else
c                  AX = AX ...
                   if(nbenrj.ge.1) then
                     MELVA1 = MCHAM1.IELVAL(1)
                     do ino1=1,nbnn1
                       MLREEL = MELVA1.IELCHE(ino1,iem1)
                       inono=inono+1
                       if(MLREEL.ne.0) then
                         RE(1,inono,NELRIG)=-1.d0*SHPP1(1,ino1)
                         RE(inono,1,NELRIG)=-1.d0*SHPP1(1,ino1)
                       endif
                     enddo
                   endif
c                  ... + B1X + B2X
                   if(nbenrj.ge.2) then
                     do jenrj=2,nbenrj
c                      on ecrit B1X pour les nbnn1 noeuds, puis B2X ...
                       MELVA1 = MCHAM1.IELVAL(jenrj)
                       do ino1=1,nbnn1
                         MLREEL = MELVA1.IELCHE(ino1,iem1)
                         if(MLREEL.eq.0) then
c                          pas d'enrichissement => on met 0
                           RX05   = 0.d0
                         else
                           SEGACT,MLREEL
                           PSIX = 0.d0
                           do iii0 = 1,nbnn1
                             PSIX = PSIX + (SHPP1(1,iii0) * PROG(iii0))
                           enddo
                           SEGDES,MLREEL
                           RX05= -1.d0*SQRT(ABS(PSIX))
                         endif
                         inono=inono+1
                         RE(1,inono,NELRIG)=RX05*SHPP1(1,ino1)
                         RE(inono,1,NELRIG)=RX05*SHPP1(1,ino1)
                       enddo
                     enddo
                   endif
                 endif

 6001          continue
 6000          continue
c           ---fin de Boucle sur les XMATRI et MELEME
c              des Partie standard et enrichie et sur les dimensions ---

c              on a trouvé le iem1 élément et on a fait le travail :
c              on passe au point de gauss suivant
               goto 1132

 1131       CONTINUE
c---------- fin de la Boucle sur les elements de structure -------------
            if(iimpi.ge.2) write(ioimp,*)
     1      'Fin de la boucle sur les elements de structure'

            if(IDEJVU(inode2).ne.10) then
         write(ioimp,*) 'Attention le noeud ',inode2,' est hors support'
            endif

 1132    CONTINUE
c======= fin de la Boucle sur les noeuds de fissure ===================

 1100 CONTINUE
c==== fin de la Boucle sur les elements de fissure =====================

c MESSAGE : Nombre de points accrochés %i1  sur  %i2 proposés
      INTERR(1)=iaccro
      INTERR(2)=NODES
      CALL ERREUR(-319)

c------------------------------------
c     AJUSTEMENT AVANT DE QUITTER
c------------------------------------

      if(iimpi.ge.2) write(ioimp,*) 'AJUSTEMENT AVANT DE QUITTER'
c     NBPTS est deja le vrai nomber de noeuds
      SEGADJ,MCOORD

C     BOUCLE SUR LES SOUS RIGIDITÉS
      ityok = 0
      DO 2000 ityty=1,(idim*NTYMAT)
         MELEME = IRIGEL(1,ityty)
         DESCR  = IRIGEL(3,ityty)
         XMATRI = IRIGEL(4,ityty)
         NBELEM = NUM(/2)
         if(NBELEM.ne.0) then
           ityok = ityok + 1
           IRIGEL(1,ityok) = MELEME
           IRIGEL(3,ityok) = DESCR
           IRIGEL(4,ityok) = XMATRI
           segdes,MELEME,XMATRI
         else
           segsup,MELEME,DESCR,XMATRI
         endif
 2000 continue
      NRIGEL = ityok
      segadj,MRIGID

c-------------------------------
c     MENAGE AVANT DE QUITTER
c-------------------------------
      segsup,MTRAV
      segdes IPT1,IPT2
      if(MCHAM1.ne.0) segdes,MCHAM1
      SEGDES,MRIGID

      if(iimpi.ge.3) write(ioimp,*) 'ecriture du MRIGID',MRIGID
      CALL ECROBJ('RIGIDITE',MRIGID)

      END