C EXCIT1    SOURCE    PV        13/09/26    21:16:37     7827           

      SUBROUTINE EXCIT1(MRIGID,MCHPO1,MCHPO2,MCHPO3,MLENTI,indic,RI2,
     &                  ipt8,ITYP)

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

*  modif : declarer inactifs les blocages redondants s'ils sont dans 
*  ipt8. Si celui-ci est non vide, on enleve les blocages inclus 
*  dedans et on s'en tient la.

-INC CCOPTIO

-INC SMRIGID
-INC SMCHPOI
-INC SMELEME
-INC SMCOORD
-INC SMLENTI

      SEGMENT SNOMIN 
        CHARACTER*4 NOMIN(0) 
      ENDSEGMENT 
      SEGMENT ISOU(IRIGEL(/2)) 
      SEGMENT ICPR(XCOOR(/1)/idimp1)
      SEGMENT JCPR(XCOOR(/1)/idimp1)
      SEGMENT KCPR(XCOOR(/1)/idimp1)
      SEGMENT YCPR(XCOOR(/1)/idimp1)
      SEGMENT ITRAV
        real*8 DVAL(NIN,NPO),DIN(NPO) 
        integer IEXI(NIN,NPO)
      endsegment
      SEGMENT  IPASS(ITEMP) 
      segment trav1
        integer iav(nbo)
        real*8 viol(nbo)
      endsegment
C
C  ***  ICPR REFERENCIE DVAL 
C  ***  NOMIN CONTIENT LES INCONNUES REFERENCES PAR LES RELATIONS 
C  ***  DVAL CONTIENT LE RESULTAT DES DEPLACEMENTS ET DES LAMBDA 
C  ***  DIN CONTIENT LES VALEURS INITIALES DU SECOND MEMBRE 
C
      logical indic
C
      idimp1 = idim+1
      segact mcoord
C 
      xcrot  = 0.
      XMCRIT = 0.
      YMCRIT = 0.
      ymult  = 1.
***   if (ityp.eq.3) ymult=1.D3
      if (ityp.eq.2) xcrot=0.250

      icpr=0
      jcpr=0
      ycpr=0
      kcpr=0
      itrav=0

* on a un maillage de relation(s) a enlever
      if (ipt8.ne.0) then
        segact,ipt8
        if (ipt8.itypel.ne.1) then
          segdes,ipt8
          call erreur(16)
          return
        endif
        if (ipt8.num(/2).ne.0) then
          segini kcpr
          do i=1,ipt8.num(/2)
            kcpr(ipt8.num(1,i))=1
          enddo
        endif
        segdes,ipt8
      endif
*
*  mchpo3 est le champ des LX sorti de excfro jcpr donne existence
*   et ycpr donne valeur
*     write (ioimp,*) ' mchpo3 dans excit1 ',mchpo3
      if (mchpo3.ne.0) then
        segini jcpr,ycpr
        segact mchpo3
        do 100 isoupo=1,mchpo3.ipchp(/1)
          msoupo=mchpo3.ipchp(isoupo)
          segact msoupo
          do 110 ic=1,nocomp(/2)
            if (nocomp(ic).eq.'LX  ') goto 120
 110      continue
          goto 140
 120      continue
          ipt2=igeoc
          mpoval=ipoval
          segact ipt2,mpoval
          do 130 iel=1,ipt2.num(/2)
            impt=ipt2.num(1,iel)
            jcpr(impt)=1
            ycpr(impt)=ycpr(impt)+vpocha(iel,ic)
 130      continue
          segdes,ipt2,mpoval
 140      continue
          segdes,msoupo
 100    continue
        segdes mchpo3
      endif
C
C** INITIALISATION DE NOMIN ET ICPR icpr repere les points appartenant 
C** a la matrice de blocage mini maxi ou ? 
C
      SEGINI,SNOMIN,ICPR
      NOMIN(**)='LX  ' 
      NPO=0 
      NBO=0 
*
      SEGACT MRIGID 
      DO 1 I=1,IRIGEL(/2) 
        ITY=IRIGEL(6,I)
        IF (ITY.EQ.0) GO TO 1
*  cas du frottement : petite verification
        if (ity.eq.2 .and. mchpo3.eq.0) then
          call erreur(721)
          return
        endif
        MELEME=IRIGEL(1,I) 
        SEGACT MELEME 
        NBO=NBO+NUM(/2) 
        IF(IIMPI.EQ.528) WRITE(IOIMP,3765) NBO 
 3765   FORMAT('  NBO  ',I5) 
        DO 2 J=1,NUM(/2) 
          DO 2 K=1,NUM(/1) 
            impt=NUM(K,J)
            IF (ICPR(impt).EQ.0) THEN
              NPO=NPO+1 
              ICPR(impt)=NPO
            ENDIF
    2   CONTINUE 
        SEGDES MELEME 
        DESCR=IRIGEL(3,I) 
        SEGACT DESCR 
        DO 3 J=1,LISINC(/2) 
          DO 4 K=1,NOMIN(/2) 
            IF (NOMIN(K).EQ.LISINC(J)) GO TO 3 
    4     CONTINUE 
          NOMIN(**)=LISINC(J) 
    3   CONTINUE 
        SEGDES DESCR 
 1    CONTINUE 
      NIN=NOMIN(/2) 
*
      IF(IIMPI.EQ.528) THEN
        WRITE(ioimp,90) NPO
        WRITE(IOIMP,9876)( NOMIN(i),i=1,NIN) 
  90    FORMAT(' ON VIENT DE PASSER LA BOUCLE 1  NPO  ',I5)
 9876   FORMAT('  NOMIN ' ,10(A4,1X)) 
      ENDIF
C 
C ****  ON REMPLIT DVAL AVEC LES VALEURS DE MCHPO2 
C   mchpo2 est le champ de deplacement propose par RESOU
C   dval(j,icpr(i)) est la j eme composante du deplacement du point i 
C   et iexi(j,icpr(i))=1
C
      SEGINI ITRAV 
c
      MCHPOI=MCHPO2 
      SEGACT MCHPOI 
      DO 5 I=1,IPCHP(/1) 
        MSOUPO=IPCHP(I) 
        SEGACT MSOUPO 
        MELEME=IGEOC 
        SEGACT MELEME 
        ITEMP=NOCOMP(/2) 
        SEGINI IPASS 
        DO 6 K=1,ITEMP
          DO 7 J=1,NIN 
            IF (NOMIN(J).EQ.NOCOMP(K)) THEN
              IPASS(K)=J 
              GO TO 6 
            ENDIF
   7      CONTINUE 
   6    CONTINUE 
        IF (IIMPI.EQ.528) WRITE(IOIMP,1555) (IPASS(KHU),KHU=1,ITEMP) 
 1555 FORMAT('   IPASS   ' ,9I10) 
        MPOVAL=IPOVAL 
        SEGACT MPOVAL 
        III=0 
        DO 8 J=1,ITEMP
          K=IPASS(J) 
          IF (K.EQ.0) GO TO 8 
          DO 9 L=1,NUM(/2) 
            IP=ICPR(NUM(1,L)) 
            IF (IP.EQ.0) GO TO 9 
*     if (k.eq.1.and.vpocha(l,j).eq.0.d0 ) then
*      write (6,*) 'LX nul dans excit1'
*      goto 9
*     endif
            DVAL(K,IP)=dval(k,ip)+VPOCHA(L,J)
            IEXI(K,IP)=1 
            III=III+1 
   9      CONTINUE 
   8    CONTINUE 
        SEGSUP IPASS 
        SEGDES MPOVAL,MSOUPO,MELEME 
   5  CONTINUE 
      SEGDES MCHPOI 
*
      IF (IIMPI.EQ.528) then
        WRITE(IOIMP,1556)III,(DVAL(1,i),i=1,NPO) 
        WRITE(IOIMP,101) 
 1556   FORMAT(' III  DVAL ',I6,/,(1X,10E12.5)) 
  101   FORMAT(' ON VIENT DE PASSER LA BOUCLE 5') 
      endif
C 
C  **** ON REMPLIT DIN PAR LES VALEURS DE MCHPO1 POUR LES LAMBDAS 
C    mchpo1 est le vecteur initial FLX ( deplacement initial)
C    din (icpr(i)) est le deplacement (FLX) initial du point i 
C
      MCHPOI=MCHPO1 
      SEGACT MCHPOI 
      DO 10 i=1,IPCHP(/1) 
        MSOUPO=IPCHP(I) 
        SEGACT MSOUPO 
        MELEME=IGEOC 
        SEGACT MELEME 
        MPOVAL=IPOVAL 
        SEGACT MPOVAL 
        DO 11 J=1,NOCOMP(/2) 
          IF (NOCOMP(J).NE.'FLX ') GO TO 11 
          DO 12 K=1,NUM(/2) 
            IP=ICPR(NUM(1,K)) 
            IF (IP.EQ.0) GO TO 12 
            DIN(IP)=VPOCHA(K,J) 
   12     CONTINUE 
   11   CONTINUE 
        SEGDES MPOVAL,MELEME,MSOUPO 
   10 CONTINUE 
      SEGDES MCHPOI 
*
      IF (IIMPI.EQ.528) then
        WRITE(IOIMP,102) 
        WRITE(IOIMP,666) (DIN(i),i=1,NPO) 
  102   FORMAT(' ON VIENT DE PASSER LA BOUCLE 10') 
  666   FORMAT('   DIN ',/,(1X,10E12.5)) 
      ENDIF
C 
C **** ON BOUCLE SUR LES RIGIDITE ET ON TESTE: 
C **** SI LE MULTI EXISTE DANS DVAL SON SIGNE PAR RAPPORT A IRIGEL(6,I) 
C **** SINON ON TESTE LA RELATION ELLE MEME A L'AIDE DES COEFF 
C **** DE LA MATRICE ET DE LA VALEUR DU LAMBDA INI (DANS DIN) 
C **** ON CREE UN TABLEAU CONTENANT LE NUMERO DES MATRICES A GARDER ET 
C **** LE NUMERO DE LA SOUS RIGIDITE 
      IPA=0 
      JG=NBO 
      SEGINI MLENTI,trav1 
      DO 313 I=1,IRIGEL(/2)
        ITY=IRIGEL(6,I)
        IF (ITY.EQ.0) GO TO 313
        MELEME=IRIGEL(1,I) 
        SEGACT MELEME 
        DESCR=IRIGEL(3,I) 
        SEGACT DESCR 
        DO 314 J=1,LISINC(/2)
          IF(LISINC(J).EQ.'LX  ') GO TO 315
  314   CONTINUE
        CALL ERREUR(5) 
        RETURN 
  315   CONTINUE
        JJ=NOELEP(J) 
        ITEMP=LISINC(/2) 
        SEGINI IPASS 
        DO 316 J=1,ITEMP
          DO 317 K=1,NIN
            IF(LISINC(J).NE.NOMIN(K)) GO TO 317
            IPASS(J)=K 
            GO TO 316
  317     CONTINUE
          CALL ERREUR (5) 
          RETURN 
  316   CONTINUE
*  RECHERCHE DU MAX DES LAMDAS POUR LE CRITERE DE DECOLLEMENT D'APPUI 
*  RECHERCHE DU MAX DES DEPLACEMENT POUR LE CRITERE DE TRAVERSEE D'APPUI
* ymcrit est le deplacement maxi de tous les points en relation unilateral
* xmcrit est le maxi des multiplicateurs existant dans le chpoin de depla
* remarque : ce chpoint de depla contient les multiplicateurs de contact
* (pression) 
        DO 30 J=1,NUM(/2) 
          IP=ICPR(NUM(JJ,J))
          DO 31 K=1,ITEMP
            IF(LISINC(K).EQ.'LX  ') GO TO 31 
            IPPP=NUM(NOELEP(K),J) 
            IPP=ICPR(IPPP) 
            YMCRIT=MAX(YMCRIT,ABS(DVAL(IPASS(K),IPP)))
   31     CONTINUE 
*  on rajoute le jeu dans le critere
          YMCRIT=MAX(YMCRIT,ABS(DIN(IP)))
          IF (IEXI(1,IP).EQ.0) GOTO 30 
          XMCRIT=MAX(XMCRIT,ABS(DVAL(1,IP))) 
  30    CONTINUE 
        SEGSUP IPASS
* on rajoute dans ymcrit la dimension de l'element
        do 32 j=1,num(/2)
          xr=0.
          yr=0.
          zr=0.
          xp=0.
          yp=0.
          zp=0.
          DO 190 K=1,ITEMP
            IF (LISINC(K).EQ.'LX  ') GO TO 190
            IPPP=NUM(NOELEP(K),J)
            ipp = idimp1*(ippp-1)
            xr = xcoor(ipp+1)
            yr = xcoor(ipp+2)
            if (idim.ge.3) zr = xcoor(ipp+3)
            goto 191
 190      continue
 191      continue
          DO 192 K=1,LISINC(/2)
            IF (LISINC(K).EQ.'LX  ') GO TO 192
            IPPP=NUM(NOELEP(K),J)
            ipp = idimp1*(ippp-1)
            xp = xcoor(ipp+1)
            yp = xcoor(ipp+2)
            if (idim.ge.3) zp=xcoor(ipp+3)
            ymcrit=max(ymcrit,abs(xp-xr),abs(yp-yr),abs(zp-zr))
 192      continue
  32    continue
 313  CONTINUE
*
      if (mchpo3.ne.0) then
*        write (6,*) ' excit1 xmcrit avant ',xmcrit
        do ip=1,ycpr(/1)
          xmcrit=max(xmcrit,abs(ycpr(ip)))
        enddo
*        write (6,*) ' excit1 xmcrit apres ',xmcrit
      endif

* test pv
*     do 320 j2=dval(/2)
*     do 320 j1=2,dval(/1)
*       ymcrit=max(abs(dval(j1,j2)),ymcrit)
* 320  continue

      XMCRIT=1D-16*XMCRIT
      YMCRIT=1D-10*YMCRIT
      ymcrit=max(xmcrit,ymcrit)
      xmcrit=xmcrit*1d-4
* Critere en cas de frottement
      yfcrit = YMCRIT*ymult
*
      DO 13 I=1,IRIGEL(/2) 
        ITY=IRIGEL(6,I)
        IF (ITY.EQ.0) GO TO 13
        MELEME=IRIGEL(1,I) 
        SEGACT MELEME 
        DESCR=IRIGEL(3,I) 
        SEGACT DESCR 
        xMATRI=IRIGEL(4,I) 
        SEGACT xMATRI 
        ITEMP=LISINC(/2) 
        SEGINI IPASS 
        DO 14 J=1,ITEMP
          IF (LISINC(J).EQ.'LX  ') GO TO 15 
   14   CONTINUE 
        CALL ERREUR(5) 
        RETURN 
   15   CONTINUE 
        JJ=NOELEP(J) 
        DO 16 J=1,ITEMP
          DO 17 K=1,NIN 
            IF(LISINC(J).NE.NOMIN(K)) GO TO 17 
            IPASS(J)=K 
            GO TO 16 
   17     CONTINUE 
          CALL ERREUR (5) 
          RETURN 
   16   CONTINUE 
        DO 18 J=1,NUM(/2) 
          IPA=IPA+1 
          IPT=NUM(JJ,J)
          IP=ICPR(IPT)
          ityz=0
*  eliminer condition redondantes pour le pas d'apres
          if (kcpr.ne.0) then
            if (kcpr(num(noelep(1),j)).ne.0) then
              goto 18
            endif
          endif
C 
C  ** CAS OU LE BLOQUAGE N'ETAIT PAS SOLLICITE 
C 
          IF (IEXI(1,IP).EQ.0) THEN 
            SS=0.D0
            DO 19 K=1,ITEMP
              if (ipass(k).eq.0) goto 19
              IF (LISINC(K).EQ.'LX  ') GO TO 19 
              IPPP=NUM(NOELEP(K),J) 
              IPP=ICPR(IPPP) 
              if (ipp.eq.0) goto 19
              SS = DVAL(IPASS(K),IPP) * RE(1,K,j) + SS
   19       CONTINUE 
            IF (IIMPI.EQ.528) WRITE(IOIMP,529) IPPP,SS 
529   FORMAT( ' LIBRE  ',I5,2X,E12.5) 
            r_z  = ABS(din(ip))*1.D-10
            r_p1 = DIN(IP) + r_z
            r_m1 = DIN(IP) - r_z
*  cas du frottement
            if (ity.eq.2) then
              ityz=jcpr(ipt)
*              write(ioimp,*) 'ipt jcpr ycpr ',ipt,jcpr(ipt),ycpr(ipt)
              if (ityz.ne.0) ityz=sign(1.1D0,ycpr(ipt))
*              if (ityz.eq.0) write (6,*) ' frottement -1 tyz ',ityz
*              write(ioimp,*) '1 dans excite ityz ',ityz
*  apparamment il faut etre plus laxiste pour le frottement
*  peut etre pas finalement
              if (ityz.eq.0) goto 20
              IF (ITYz.EQ. 1.AND.SS.LE.r_p1+yfcrit) GOTO 20
              IF (ITYz.EQ.-1.AND.SS.GE.r_m1-yfcrit) GOTO 20
            endif
            IF (ITY.EQ. 1.AND.SS.LE.r_p1+YMCRIT) GOTO 20
            IF (ITY.EQ.-1.AND.SS.GE.r_m1-YMCRIT) GOTO 20
            LECT(IPA)=1 
            viol(ipa)=ss-din(ip)-ymcrit
       if (iimpi.eq.1967)
     >     write (6,*) ' ss din ymcrit nouveau blocage '
     $   ,ss,din(ip),ymcrit,i,ipa,viol(ipa),ipt,ityz,ityp
*     on a un (1) nouveau blocage on arrete
   20       CONTINUE 
C 
C  ** CAS OU LE BLOQUAGE ETAIT SOLLICITE  PETIT PROBLEME DE TEST DE 
C     PRECISION SUR LA VALEUR DE LA REACTION 
C 
          ELSE 
            iav(ipa)=1
            SS=DVAL(1,IP) 
            IF(IIMPI.EQ.528) WRITE(IOIMP,530) NUM(JJ,J),SS 
  530 FORMAT('  BLOQUER ' ,I5,2X,E12.5) 
*  cas du frottement
            if (ity.eq.2) then
              ityz=jcpr(ipt)
              if (ityz.ne.0) ityz=sign(1.1D0,ycpr(ipt))
*              if (ityz.eq.0) write (6,*) ' frottement -2 ityz ',ityz
*  apparamment il faut etre plus laxiste pour le frottement
*  peut etre pas finalement
              if (ityz.eq.0) goto 21
              IF(ITYz.EQ.1.AND.SS.LT. XMCRIT) GO TO 21
              IF(ITYz.EQ.-1.AND.SS.GT.-XMCRIT) GO TO 21
            endif
            IF(ITY.EQ.1.AND.SS.LT. XMCRIT) GO TO 21
            IF(ITY.EQ.-1.AND.SS.GT.-XMCRIT) GO TO 21
            LECT(IPA)=1 
            goto 18
   21       CONTINUE 
            viol(ipa)=ss+xmcrit
C      write (6,*) ' ss xmcrit ',ss,xmcrit,ityz,i,ipa,
C     >  num(3,j),num(4,j),num(5,j)
      if (iimpi.eq.1967)
     > write (6,*) ' appui disparait ',ipa,ss,xmcrit,ipt,ityz,ityp
          ENDIF 
   18   CONTINUE 
        SEGSUP IPASS 
        SEGDES MELEME,DESCR,xMATRI 
   13 CONTINUE 
      IF(IIMPI.EQ.528) WRITE(IOIMP,*) 'ON VIENT DE PASSER LA BOUCLE 13'
C 
C  **** DANS LECT ON DIT SI LA JEEME RIGI ELEMTAIRE EST A CONSERVER 
C 
      NRIGEL=IRIGEL(/2) 
      SEGINI RI2
C 
C  **** CAS OU IL N'Y A RIEN A GARDER  ON CREE UNE RIGIDITE VIDE 
C  **** POUR POUVOIR CONTINUER D'ITERER SI IL Y A LIEU 
C 
      IF (NRIGEL.EQ.0) THEN
        IF (IIMPI.EQ.528) WRITE(IOIMP,*) ' IL N''Y A RIEN A CREER'
        GO TO 50
      ENDIF 
*
*** ne garder que ce qui viole de plus de xcrot du max
*
* recherche du max
      violmf=0.
      violmd=0.
      imf=0
      imd=0
      do 40 i=1,nbo
        violm=abs(viol(i))
        if (iav(i).eq.0.and.lect(i).ne.0) then
          if (violm.gt.violmd) then
            violmd=violm
            imd=i
          endif
        else if (iav(i).eq.1.and.lect(i).ne.1) then
          if (violm.gt.violmf) then
            violmf=violm
            imf=i
          endif
        endif
  40  continue
*      write (6,*) ' ityp xcrot violmd violmf ',ityp,xcrot,
*     >            violmd,violmf,imd,imf
****  if (ityp.eq.3) then
****    do 42 i=1,nbo
****      if (iav(i).eq.0.and.imd.ne.i) lect(i)=iav(i)
****      if (iav(i).eq.1.and.(imd.ne.0.or.imf.ne.i)) lect(i)=iav(i)
*42     continue
****  else
**      rvd = xcrot*violmd
**      rvf = xcrot*violmf
        rvd = 0.250*violmd
        rvf = 0.707*violmf
        if (ityp.eq.3) then
         rvd=rvd*1d-14
         rvf=rvf*1d-14
        endif
        do 41 i=1,nbo
          violm=abs(viol(i))
          if (imf.ne.0.and.iav(i).eq.0) lect(i)=iav(i)
**        if (imd.ne.0.and.iav(i).eq.1) lect(i)=iav(i)
          if (iav(i).eq.0.and.violm.ge.rvd) goto 41
**?          if (iav(i).eq.1.and.violm.ge.rvf) goto 41
          if (iav(i).eq.1) goto 41
          lect(i)=iav(i)
  41    continue
****  endif
C 
C  ** IL FAUT CREER UNE RIGIDITE 
C 
      IPA= 0 
      IRI=0 
      DO 25 I =1,IRIGEL(/2) 
        IF(IRIGEL(6,I).EQ.0) GO TO 25
        MELEME=IRIGEL(1,I) 
        SEGACT MELEME 
        NELRIG=0 
        DO 27 J=1,NUM(/2) 
          IF(LECT(IPA+J).EQ.1) NELRIG=NELRIG+1 
   27   CONTINUE 
        if (nelrig.eq.0) goto 26
        IRI=IRI+1
        xMATRI=IRIGEL(4,I) 
        SEGACT xMATRI 
        nligrp=re(/2)
        nligrd=re(/1) 
        SEGINI xMATR1 
        RI2.IRIGEL(4,IRI)=xMATR1 
        DESCR=IRIGEL(3,I) 
        SEGINI,DES1=DESCR 
        SEGDES DES1 
        RI2.IRIGEL(3,IRI)=DES1 
        RI2.IRIGEL(5,IRI)=IRIGEL(5,I) 
        RI2.IRIGEL(6,IRI)=IRIGEL(6,I) 
        RI2.COERIG(IRI)=COERIG(I) 
        NBNN=NUM(/1) 
        NBELEM=NELRIG 
        NBSOUS=0 
        NBREF=0 
        SEGINI IPT1 
        IPT1.ITYPEL=ITYPEL 
        RI2.IRIGEL(1,IRI)=IPT1 
        I2=0 
        DO 28 J=1,NUM(/2) 
          IF(LECT(IPA+J).EQ.0) GO TO 28 
          I2=I2+1 
          DO 29 K=1,NUM(/1) 
            IPT1.NUM(K,I2)=NUM(K,J) 
 29       CONTINUE 
          do io=1,nligrp
            do iu=1,nligrd
              xmatr1.re(iu,io,i2)=re(iu,io,j)
            enddo
          enddo 
 28     CONTINUE 
        SEGDES xMATR1,xMATRI,IPT1 
 26     CONTINUE 
        IPA=IPA+NUM(/2) 
        SEGDES MELEME 
   25 CONTINUE 
*
      if (iri.ne.ri2.irigel(/2)) then
         nrigel=iri
         segadj ri2
      endif

   50 CONTINUE 
      SEGDES,RI2,MRIGID 

*  indice de retour
      indic = .true.
      do 55 i = 1, nbo
       if (iav(i).ne.lect(i)) indic = .false.
  55  continue
      SEGDES MLENTI

      SEGSUP SNOMIN,ICPR,ITRAV,trav1
      if (mchpo3.ne.0) segsup jcpr,ycpr
      if (kcpr.ne.0) segsup kcpr 

      RETURN 
      END