C ACCRO     SOURCE    BP208322  12/12/17    21:15:06     7615           
      SUBROUTINE ACCRO(igliss)
C========================================================================      
C   ACCROCHAGE D UN MAILLAGE  SUR UN MASSIF PAR L INTERMEDIAIRE 
C   DES FONCTIONS DE FORME 
C   SUBROUTINE APPELLEE PAR RELA2                 
C   JMB      OCTOBRE  1996    
C========================================================================  
                                                               
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)  
-INC CCREEL                                                
-INC CCOPTIO                                                                    
-INC SMCOORD                                                                    
-INC SMELEME 
-INC SMRIGID
-INC SMLMOTS
-INC CCHAMP
-INC CCGEOME
-INC SMMODEL
      external shape
      DATA EPSI/1.D-9/
      DIMENSION ICOR(100),IMEL(100),imtt(100)
      SEGMENT IEXCLU
      INTEGER IDEJVU(NBPTI)       
      ENDSEGMENT             
      SEGMENT ICOUNT 
      INTEGER IEINT(2,NODES),MAXPZ(1)
      REAL*8 QQQ(3,NODES)
      ENDSEGMENT
      SEGMENT ICOSOM
      INTEGER IPTT(NBSZ)
      ENDSEGMENT   
      SEGMENT MTRAV 
      REAL*8 CORDEL(3,NBNN),QSI(3),XPO(3),SHPP(6,NBNN),COE(IDIM,30)
      REAL*8 CC(100),XEL(3,NBNN),BPSS(3,3),XPT(3),CORED(3,NBNN)
      ENDSEGMENT
      SEGMENT IOUBLI
      INTEGER IPTOU(NBOUB,2)
      ENDSEGMENT
      segment itange
      real*8 xperp1(idim,nbpti),xperp2(idim,nbpti)
      integer iextr(nbpti)
      endsegment      
      segment inomil(0) 
      segment icpr(xcoor(/1)/(idim+1))
      segment imul(idim,nbpt)
      segment icomu(nbpt)
      logical ltelq


C LECTURE des MAILLAGE et MMODEL ***********************************     

      IFLAG= 0                                                            
      IFROCA=0
      MELEME=0                                                       
      CALL LIROBJ('MAILLAGE',IPP2,1,IRETOU)
      IF(IERR.NE.0) RETURN 
      CALL LIROBJ('MAILLAGE',MELEME,0,IRETOU)
      MELEMI=MELEME
      IF(IRETOU.EQ.0) then       
         melemi=ipp2
         meleme=ipp2     
         CALL LIROBJ('MMODEL',IPMODL,1,IFROCA)                          
         IF (IERR.NE.0) RETURN  
      endif    
c     ! Erwan


C**** CAS OU UN MODELE A ETE LU **************************************
      if (IFROCA.ne.0) then                                            
        mmode1 = ipmodl
        segact mmode1
        nbpt=0
        segini icpr
C creation des l'ensemble des points a accrocher et des multiplicateurs
C correspondant. imul(i,j) est le ieme multiplicateur a accrocher au j eme 
C point local ( icpr(ia)=j) 
C --- puis Assemblage du maillage global (tous les câbles)
        do 500 isouf = 1,mmode1.kmodel(/1)
          imode1 = mmode1.kmodel(isouf)
          segact imode1
          if(imode1.tymode(1).eq.'MMODEL') then
            mmode2 = imode1.ivamod(1)
            segact mmode2
            imode2 = mmode2.kmodel(1)
          else
            call erreur(975)
            return
          endif                 
          segact imode2
          ipt7 = imode2.imamod
          if(isouf.eq.1) then
            ipt5=ipt7
          else
            ltelq=.false.
            call fuse(ipp4,ipt7,ipt5,ltelq)
            ipt5=ipp4
          endif
          ipt6=imode1.imamod
          segact ipt6
          if(ipt6.itypel.ne.22) then
             call erreur (19)
             return
          endif 
          segdes imode2,mmode2
          do 501 iel=1,ipt6.num(/2)
             ia=ipt6.num(2,iel)
             if(icpr(ia).eq.0 ) then
              nbpt=nbpt+1
              icpr(ia)=nbpt
             endif
  501     continue
  500   continue         
        segini imul,icomu
        nbpt=0        
        do 502 isouf = 1,mmode1.kmodel(/1)
          imode1 = mmode1.kmodel(isouf)
          ipt6=imode1.imamod
          do 503 iel=1,ipt6.num(/2)
             ia=ipt6.num(2,iel)
             ib=icpr(ia)
             icomu(ib)=icomu(ib)+1                       
             imul(icomu(ib),ib)=ipt6.num(1,iel)
  503     continue
          segdes imode1,ipt6
  502   continue 
        segdes mmode1          
        ipp2=ipt5
*  fin du if(ifroca..;
      endif
C**** FIN DU PRE-TRAITEMENT DANS LE CAS D'UN MODELE ******************


C**** LECTURE DES NOMS D'INCONNUE ET DU CRITERE *********************

      CALL LIROBJ('LISTMOTS',MLMOTX,0,IRETOU)
      CCRIT =1.D-5 
      CALL LIRREE(CCRIT,0,IRETOU)
      CCRIT = -ABS(CCRIT) 
      
      IF(MLMOTX.EQ.0) THEN  
C     options par defaut
        JGN=4
        JGM=IDIM
        NINC= IDIM        
        SEGINI MLMOTS 
        IF(IFOUR.EQ.0)  THEN
C       axisymetrique 
           MOTS(1)='UR  '
           MOTS(2)='UZ  '
        ELSE         
           MOTS(1)='UX  '
           MOTS(2)='UY  '
           IF (IDIM.EQ.3) MOTS(3)='UZ  '
        ENDIF
      ELSE
C     on a lu les noms d'inconnues
           if(igliss.eq.1) then
             call erreur (19)
             return
           endif 
           MLMOTS = MLMOTX 
           SEGACT MLMOTS 
           NINC = MOTS(/2)                      
           IF(NINC.GT.100) THEN
             SEGDES MLMOTS
             CALL ERREUR(19)
             RETURN
           ENDIF
      ENDIF 
C     IF (MLMOTX.EQ.0) THEN
C      MOTERR(1:8)='LISTMOTS'
C      CALL ERREUR(37)
C       RETURN
C     ENDIF 


C**** VERIFICATION DU MAILLAGE MELEME ********************************

C   LE SECOND  MAILLAGE DOIT ETRE MASSIF 
C   POUR LE MOMENT  TETRAEDRES ET PYRAMIDES SONT EXCLUS 
      SEGACT MELEME 
      IPT1 = MELEME
      DO 19 IOB=1,MAX(1,LISOUS(/1))                                             
        IF (LISOUS(/1).NE.0) THEN                                                 
          IPT1=LISOUS(IOB)
          SEGACT IPT1 
        ENDIF                                                          
        ITY = IPT1.ITYPEL 
        IF ((idim.eq.2.and.ITY.LT.4).or.
     &      (idim.eq.3.and.ity.lt.14)) THEN
          CALL ERREUR(16)
          SEGDES IPT1*NOMOD         
          RETURN    
        ENDIF
        SEGDES IPT1*NOMOD                                                                  
19    CONTINUE    


C**** TRAVAIL ********************************************************

      SEGACT MCOORD 
      NBPTI = XCOOR(/1)/(IDIM+1)
C on veut savoir si il y a des point vus deux fois (interfaces sous zones)       
      SEGINI IEXCLU
      IEXX = IEXCLU 
           
      IPT2 = IPP2 
      SEGACT IPT2
C  ON VA CHANGER LE MAILLAGE LINEAIRE EN ELEMENTS POI1 
C   SI IL NE N  L EST PAS DEJA   
      IF((IPT2.ITYPEL).NE.1) THEN            
        CALL CHANGE(IPP2,1) 
      ENDIF 
      SEGDES IPT2
      
C--------------------------------------        
*---- cas de l'accrochage glissant ----
      IF ( igliss.eq.1) THEN
      
*       on va attacher à chaque point un vecteur tangent au cable
*       si on est arrive avec un modele, on prend le maillage du modele 
*       en reference
        segini itange
        idim1=idim+1
        segact ipt2
        meleme = ipt2
c       boucle sur les maillages elementaires        
        do 3657 i=1,max(1,ipt2.lisous(/1))
          if(ipt2.lisous(/1).ne.0) then
            meleme= ipt2.lisous(i)
            segact meleme
          endif
          if(itypel.ne.2) then
            call erreur (946)
            return
          endif
c         boucle sur les elements
          do 3658 j=1,num(/2)
            ia=num(1,j)
            ib=num(2,j)
            if(iextr(ia).gt.2.Or.iextr(ib).gt.2) then
              call erreur (947)
              return
            endif 
            xab=xcoor((ib-1)*idim1 + 1)- xcoor((ia-1)*idim1 + 1) 
            yab=xcoor((ib-1)*idim1 + 2)- xcoor((ia-1)*idim1 + 2)
            zab=0.
            if(idim.eq.3) then
              zab = xcoor((ib-1)*idim1 + 3)- xcoor((ia-1)*idim1 + 3)
            endif
            xlo=sqrt(xab * xab + yab * yab + zab * zab)
            if(xlo.eq.0.) then
             call erreur(945)
             return
            endif
            xab=xab/xlo
            yab=yab/xlo
            zab=zab/xlo  
            if(iextr(ia).eq.0) then
              xperp1(1,ia)=xab
              xperp1(2,ia)=yab
              if(idim.eq.3) xperp1(3,ia)=zab
            else 
              xxx=1.d0
              if( iextr(ia).eq.1) xxx = -1.d0
              xperp1(1,ia)=xab*xxx+ xperp1(1,ia)
              xperp1(2,ia)=yab*xxx+ xperp1(2,ia)
              if(idim.eq.3) xperp1(3,ia)=zab*xxx+ xperp1(3,ia)
              xlo=sqrt (xperp1(1,ia)**2+xperp1(2,ia)**2)
              if(idim.eq.3)xlo=sqrt (xperp1(1,ia)**2+xperp1(2,ia)**2
     $        + xperp1(3,ia)**2)
              xperp1(1,ia)=xperp1(1,ia)/xlo
              xperp1(2,ia)=xperp1(2,ia)/xlo
              if(idim.eq.3)xperp1(3,ia)=xperp1(3,ia)/xlo
            endif
            if(iextr(ib).eq.0) then
              xperp1(1,ib)=xab
              xperp1(2,ib)=yab
              if(idim.eq.3) xperp1(3,ib)=zab
            else
              xxx=1.d0 
              if( iextr(ia).eq.2) xxx = -1.d0      
              xperp1(1,ib)=xab*xxx+ xperp1(1,ib)
              xperp1(2,ib)=yab*xxx+ xperp1(2,ib)
              if(idim.eq.3) xperp1(3,ib)=zab*xxx+ xperp1(3,ib)
              xlo=sqrt (xperp1(1,ib)**2+xperp1(2,ib)**2)
              if(idim.eq.3) xlo=sqrt (xperp1(1,ib)**2+xperp1(2,ib)**2
     $        + xperp1(3,ib)**2)
              xperp1(1,ib)=xperp1(1,ib)/xlo
              xperp1(2,ib)=xperp1(2,ib)/xlo
              if(idim.eq.3)xperp1(3,ib)=xperp1(3,ib)/xlo               
            endif
            iextr(ia)= iextr(ia)+1
            iextr(ib)= iextr(ib)+2
 3658     continue
c         fin de boucle sur les elements        
          if(ipt2.lisous(/1).ne.0) segdes meleme       
 3657   continue
c       fin de boucle sur les maillages elementaires        

        segdes ipt2
        do 3659 i=1,nbpti
          if(iextr(i).eq.0)go to 3659
          if(abs(xperp1(1,i)).gt.2.e-2.OR.abs(xperp1(2,i)).gt.2.e-2)
     $    then
            xperp2(1,i)=-xperp1(2,i)
            xperp2(2,i)=xperp1(1,i)
          else
            xperp2(1,i)=-xperp1(3,i)
            xperp2(3,i)= xperp1(1,i)
          endif
          if(idim.eq.3) then        
            xlo=sqrt (xperp2(1,i)**2+xperp2(2,i)**2
     $      + xperp2(3,i)**2)
            xperp2(1,i)=xperp2(1,i)/xlo
            xperp2(2,i)=xperp2(2,i)/xlo
            xperp2(3,i)=xperp2(3,i)/xlo
            xaa = xperp1(3,i)*xperp2(2,i)-xperp1(2,i)*xperp2(3,i)
            xbb = xperp1(1,i)*xperp2(3,i)-xperp1(3,i)*xperp2(1,i)
            xcc = xperp1(2,i)*xperp2(1,i)-xperp1(1,i)*xperp2(2,i)
            xperp1(1,i)=xaa
            xperp1(2,i)=xbb
            xperp1(3,i)=xcc
*          write(6,*) 'xaa xbb xcc',xaa,xbb,xcc
          else
            xperp1(1,i)=xperp2(1,i)
            xperp1(2,i)=xperp2(2,i)
          endif
 3659   continue
      ENDIF 
*---- fin du traitement specifique accrochage glissant ----
C----------------------------------------------------------        

      IPT2 = IPP2 
      SEGACT IPT2                    
      NODES = IPT2.NUM(/2)   
      SEGDES IPT2 
C 
      MELEME=MELEMI         
      SEGACT MELEME 
      IPP0 = MELEME    
      NBSZ = MAX(1,LISOUS(/1))
      IPP1 = MELEME 
      SEGINI ICOSOM
      iaccro = 0 

C---- BOUCLE SUR LES SOUS ZONES ----                                                
      DO 22 IOB=1,MAX(1,LISOUS(/1))                                             
      IF (LISOUS(/1).NE.0) THEN                                                 
       IPP1=LISOUS(IOB)                                                                                                                           
      ENDIF             
C     ICI LE  ZONAGE 
C============================================
      CALL  ZONAG(IPP1,IPP2,ICC,IEXX,0) 
C============================================
C     TABLE DES POINTEURS SUR LES TABLEAUX DE CORRESPONDANCE 
C     WRITE(6,*) 'POINTEUR ZONE ',IOB,ICC 
      IPTT(IOB) = ICC
      ICOUNT = ICC
      SEGACT ICOUNT
      iaccro = iaccro + MAXPZ(1)
c       do 779 i=1,MAXPZ(1)
C       write(6,*)' i ient   ', i,IEINT(1,I),IEINT(2,I)
c  779  continue       
      SEGDES ICOUNT
22    CONTINUE 
C---- FIN DE LA BOUCLE SUR LES SOUS ZONES ----                                                

C---- TABLEAU  des POINT PAS VUS ----
      NBOUB  = NODES - IACCRO 
C     write(6,*) 'nombre de points exterieurs ' ,nboub
      IPT2=IPP2  
      SEGACT IPT2    
      SEGINI IOUBLI
      SEGACT IEXCLU
      IUO = 0 
      DO 21 I=1,NODES 
       IP = IPT2.NUM(1,I)
       IF(IDEJVU(IP).EQ.0) THEN
         IUO = IUO + 1 
         IPTOU(IUO,1)=IP        
       ENDIF
21    CONTINUE 

C-----------------------------------------------------------         
C---- RATTRAPAGE EVENTUEL DE POINTS SUR LA LIMITE ----
      IF(IUO.NE.0) THEN
      
      MELEME=IPP0
      SEGACT MELEME 
      IPT1=MELEME 
      krat = 0 
      
C---- boucle sur les sous zones ----
      DO 1010 IOB=1,MAX(1,LISOUS(/1))
        IF (LISOUS(/1).NE.0) THEN
          IPT1=LISOUS(IOB)
          segact ipt1
        ENDIF
        
        iele=ipt1.itypel      
        kd = kdegre(iele) 
        nso = nbsom(iele)
        iad=nspos(iele)-1
        segini inomil
          if(kd.eq.3) then
C element quadratiques        
             do 762  i=1,ipt1.num(/1)     
                do 763  j=1,nso 
                   iso = ibsom(iad+j)        
                   if(i.eq.iso) goto 762 
 763            continue
                inomil(**)= i  
 762         continue    
          elseif(kd.eq.2) then
C  elements  lineaires        
             do i=1,ipt1.num(/1) 
                inomil(**)=i
             enddo        
          endif 
        NBNN = ipt1.num(/1)
        SEGINI MTRAV 
        ICOUNT= IPTT(IOB)
        SEGACT ICOUNT*MOD      
3333    FORMAT(I4,I6,4X,8I4)   
3338    format(I8,3E12.5)
C                    
C...... boucle sur les elements de la sz .......
        DO 1121 iem=1,ipt1.NUM(/2)
        call doxe(xcoor,idim,nbnn,ipt1.num,iem,xel) 
        
C.......boucle sur les points candidats au rattrapage ......
        DO 1120 IPT = 1,IUO 
        
          IP = IPTOU(IPT,1)
          if(ip.LT.0) GOTO 1120
          do jl = 1,ipt1.num(/1)
            if (ip.eq.ipt1.num(jl,iem)) goto 1120 
          enddo 
C
          IREFP= (IP-1)*(IDIM+1)+1                 
          XPO(1)= XCOOR(IREFP)
          XPO(2)= XCOOR(IREFP+1)
          XPO(3)= 0.D0
          IF (IDIM.EQ.3) XPO(3)=XCOOR(IREFP+2)            
C         coordonnes intrinseques et fonctions de forme : qsijs
C         write(6,5777 ) ((xel(ia,ib),ia=1,idim ),ib=1,nbnn)
5777      format(6e12.5)              
          call qsijs(xel,iele,nbnn,idim,XPO,SHPP,qsi,iret)     
          IF (IRET.NE.0) goto 1120
C         CALL ERREUR(797)
C         RETURN
C         ENDIF 
          do i=1,inomil(/1)
            ilp= inomil(i)
            if( shpp(1,ilp).lt.ccrit) goto 1120
          enddo
C         ce point exterieur  est rattrapable 
          Krat = krat+1 
          NNP = MAXPZ(1)+1
          MAXPZ(1)= NNP
          IEINT(1,NNP) = ip
          iptou(ipt,1) =-ip          
          IEINT(2,NNP) = iem 
          do i=1,idim
            QQQ(i,NNP) =qsi(i) 
          enddo   
C       write(6,2375)ip,iem,(XPO(i1),i1=1,idim),(qsi(i2),i2=1,idim)

1120    continue 
C.......fin de la boucle sur les points candidats au rattrapage ......
        
2375    format(2i4,6e12.5)              
1121    CONTINUE
C...... fin de la boucle sur les elements enveloppe de la sz .......
        SEGDES ICOUNT
        segsup inomil,MTRAV
        IF (LISOUS(/1).NE.0) SEGDES IPT1  
        
1010  continue
C---- fin de la boucle sur les zones ---- 
       
C     comptage du nombre de points rattrapes 
      k= 0
      do 777 i=1,IUO
        IF(IPTOU(I,1).GT.0) GOTO 777
        k= k+1         
777   continue
C
C     write(6,*) ' nombre de points rattrapes ' ,k        
C     on donne le nombre de points non projetés
C     if (k.ne.iuo) then 
C       INTERR(1)=iuo-k
C       CALL ERREUR(-322)   
C     endif
      iaccro = iaccro + k      
       
      ENDIF  
C---- FIN DE LA SEQUENCE DE RATTRAPAGE         
C-----------------------------------------------------------         

c MESSAGE : Nombre de points accrochés %i1  sur  %i2 proposés
      INTERR(1)=iaccro
      INTERR(2)=NODES
      CALL ERREUR(-319)
C
C
      MELEME = IPP0 
      SEGACT MELEME 
      IPT1= MELEME           
      NRIGEL = 0
      NBSZ = 0
      NBPTOT = 0  
      if( igliss.eq.1.and .ifroca.eq.0) ninc=ninc-1
     
C  NOMBRE DE SOUS ZONES TOUCHÉES  ET NOMBRE TOTAL DE NOEUDS CONCERNÉS 
      DO 399  IOB=1,MAX(1,LISOUS(/1))
      ICOUNT= IPTT(IOB)
C        write(6,*)  ' segment active' ,IPTT(IOB)
      SEGACT ICOUNT
      IF(MAXPZ(1).EQ.0) GOTO 398      
      IF(MAXPZ(1).EQ.0) GOTO 398
      NBSZ = NBSZ + 1
      NBPTOT = NBPTOT + MAXPZ(1) 
398   SEGDES ICOUNT                       
399   CONTINUE

  
C     dans le cas ifroca ne 0 les multiplicateurs existent deja
      if(ifroca.eq.0) then     
        NBPTS = NBPTI + NBPTOT*NINC
        SEGADJ MCOORD
      endif
      NRIGEL = NBSZ * NINC
C     INITIALISATION  DE MRIGID 
      NRIGE = 8
      SEGINI MRIGID 
      IFORIG = IFOMOD
      MTYMAT ='RIGIDITE' 
      ICZ = 0 
      INRI=0 
C________________________________________________________________
C____ BOUCLE SUR LES SOUS ZONES __________________________                                                
      DO 400 IOB=1,MAX(1,LISOUS(/1))
      
      ICOUNT = IPTT(IOB)
      SEGACT ICOUNT             
      IF (MAXPZ(1).EQ.0) GOTO 395 
      ICZ = ICZ + 1 
C      WRITE(6,*) 'SOUS ZONE ',IOB ,'NOMBRE DE POINTS ',MAXPZ(1),ICOUNT      
      IF (LISOUS(/1).NE.0) THEN                                                 
       IPT1=LISOUS(IOB)                                                         
       SEGACT IPT1                                                              
      ENDIF
      KEL = IPT1.ITYPEL 
      NBNN = IPT1.NUM(/1) 
C      
      SEGINI MTRAV 
C   ce sera ici               
      NBELEM  =  MAXPZ(1)     
      NBMA1 =  IPT1.NUM(/1)      
      NBNN = NBMA1 + 2
      NBSOUS = 0
      NBREF = 0 
C     creation des xmatri
      NELRIG  = NBELEM 
      NLIGRD  = IPT1.NUM(/1)+1+1     
      if(igliss.eq.1) nligrd = 1 + idim*(IPT1.NUM(/1)+1)
      NLIGRP =  NLIGRD 
      do iop=1,ninc
        SEGINI XMATRI 
        imtt(iop)=Xmatri
      enddo
C      
      INRI =NINC *(ICZ-1)+1
C      write(6,*) ' ninc' , ninc
c     boucle sur les xmatri                                                
      DO  201 I= 1,NINC
        SEGINI IPT3 
        IPT3.ITYPEL= 22  
        if( ifroca.NE.0.AND.i.eq.ninc)  irigel(6,INRI+I-1)=2
C      write (6,*) ' ipt3 ',ipt3.num(/1),ipt3.num(/2),ipt1.num(/1),
C     >  ipt1.num(/2)
        IMEL(I)= IPT3
        IRIGEL (1,INRI+I-1) = IPT3
        IRIGEL (4,INRI+I-1) = IMtt(i)
        COERIG(INRI+I-1)= 1.D0
201   CONTINUE       
      NELRIG = 0         
      
C________________________________________________________________ 

C.... boucle sur les elements .......
      DO 300 K = 1,NBELEM       
C       NELRIG = NELRIG + 1
        IP  = IEINT(1,K)
        IEL = IEINT(2,K)
        irefp=(ip-1)*(idim+1)+1
        xpo(1)=xcoor(irefp)
        xpo(2)=xcoor(irefp+1)
        xpo(3)= 0.D0
        if(idim.eq.3) xpo(3)=xcoor(irefp+2) 
        qsi(1)=qqq(1,k)
        qsi(2)=qqq(2,k)     
        qsi(3)=qqq(3,k)      
    
        CALL SHAPE(QSI(1),QSI(2),QSI(3),KEL,SHPP,IRET)                         
      
c      -boucle sur les xmatri-
        DO 59 IK =1,NINC
      
        ipt3=imel(ik)
C       write(6,*) ' IPT1.NUM(/1) ipt1;num(/2)ipt3.num(/1) ipt3.num(/2)'
C       write(6,*) IPT1.NUM(/1), ipt1.num(/2),ipt3.num(/1), ipt3.num(/2)  
        if(ifroca.eq.0) then    
          NBPTI = NBPTI + 1 
          IDE =( NBPTI-1)*(IDIM+1 )
          XCOOR(IDE+1) = XPO(1) 
          XCOOR(IDE+2) = XPO(2)      
          IF(IDIM.EQ.3) THEN 
           XCOOR(IDE+3) = XPO(3) 
          ENDIF 
        else
         ibp=icpr(ip)
         nbpti= imul(ik,ibp)
        endif
      
C       REMPLISSAGE DU MELEME IPT3
        IPT3.NUM(1,K) = NBPTI
        IPT3.NUM(2,K) = IP 
        DO 301 LI = 1,IPT1.NUM(/1)
          IPT3.NUM(2+LI,K) = IPT1.NUM(LI,IEL) 
301     CONTINUE

C       REMPLISSAGE DU XMATRI
        xmatri=imtt(ik) 
        CC(1)   = 1.D0
        NHD = IPT1.NUM(/1)+1 
        DO 305 KM=1,IPT1.NUM(/1) 
           CC(1+KM) = -SHPP(1,KM)  
305     CONTINUE 

c       Cas glissant                       
        IF(igliss.eq.1) THEN
         xaa = xperp1(1,ip)
         xbb = xperp1(2,ip)
         if(idim.eq.3) xcc = xperp1(3,ip)
         if( ik.eq.2) then
           xaa = xperp2(1,ip)
           xbb=xperp2(2,ip)
           if(idim.eq.3) xcc= xperp2(3,ip)
         endif
         if(ifroca.ne.0.and.ik.eq.ninc) then
C petit calcul du vercteur tangent
           if( idim.eq.2) then
             XAA = xperp1(2,ip) 
             XBB =-xperp1(1,ip)
           else
             XAA = xperp1(2,ip)*xperp2(3,ip)-xperp1(3,ip)*xperp2(2,ip) 
             XBB = xperp1(3,ip)*xperp2(1,ip)-xperp1(1,ip)*xperp2(3,ip)   
             XCC = xperp1(1,ip)*xperp2(2,ip)-xperp1(2,ip)*xperp2(1,ip)
           endif
         endif           
         DO 337 km=1,nhd
         ccC = -cc(KM)
         RE(1,1+(km-1)*IDIM +1,k)= -CCC*XAA
         RE(1,1+(km-1)*IDIM +2,k)= -CCC*XBB
         RE(1+(km-1)*IDIM +1,1,k)= -CCC*XAA
         RE(1+(km-1)*IDIM +2,1,k)= -CCC*XBB
         IF(IDIM.EQ.3) THEN
           RE(1,1+(km-1)*IDIM +3,k)= -CCC*XCC
           RE(1+(km-1)*IDIM +3,1,k)= -CCC*XCC
         ENDIF
  337    CONTINUE
  
c       Cas non glissant                       
        ELSE                  
         DO 306 KM =1,NHD
         CCC = CC(KM)
         RE(1,KM+1,k)= CCC     
         RE(1+KM,1,k)= CCC    
306      CONTINUE 
        ENDIF

  59    continue       
c     - fin de boucle sur les xmatri -

300   CONTINUE 
C.... fin de boucle sur les elements .......
      SEGDES xMATRI
      do ik=1,ninc
        IPT3 = IMEL(IK)
        SEGDES IPT3  
      enddo           

C     identification des duals  des ddl  en question
c     NL = 20 
c bp 18.01.2011 : mieux vaut utiliser LNOMDD renseigné dans cchamp+bdata
      NL = LNOMDD
      DO 326 I = 1,MOTS(/2)
      DO 327 KM = 1,NL
       IF(NOMDD(KM).EQ.MOTS(I)) THEN 
       ICOR(I) = KM 
       GOTO 326 
       ENDIF
327   CONTINUE       
326   CONTINUE   

      NLIGRD  = IPT1.NUM(/1)+1+1     
      if(igliss.eq.1) nligrd = 1 + idim*(IPT1.NUM(/1)+1)
      NLIGRP =  NLIGRD 
c---- boucle sur les xmatri ----
      DO 328 I=1,NINC 
      
       NLIGRP =  NLIGRD          
       SEGINI DESCR 
       LISINC(1)= 'LX' 
       LISDUA(1)= 'FLX'
       NOELEP(1)= 1  
       NOELED(1)= 1  
       
c    - Cas non glissant - 
       if(igliss.eq.0) then
         DO 307 KM1 = 1,IPT1.NUM(/1)+1
           LISINC(KM1+1)= MOTS(I)
           LISDUA(KM1+1)= NOMDU(ICOR(I))                             
           NOELEP(KM1+1)= KM1+1              
           NOELED(KM1+1)= KM1+1           
307      CONTINUE 

c    - Cas glissant - 
       else
         nhd=IPT1.NUM(/1)+1
         do 347 km1=1, nhd
           lisinc(1+(km1-1)*idim +1)=MOTS(1)
           lisinc(1+(km1-1)*idim +2)=MOTS(2)   
           LISDUA(1 +(km1-1)*idim +1)= NOMDU(ICOR(1))    
           LISDUA(1 +(km1-1)*idim +2)= NOMDU(ICOR(2))
           NOELEP(1+(km1-1)*idim +1)= KM1+1
           NOELEP(1+(km1-1)*idim +2)= KM1+1 
           NOELED(1+(km1-1)*idim +1)= KM1+1 
           NOELED(1+(km1-1)*idim +2)= KM1+1       
           if(idim.eq.3) then
             lisinc(1+(km1-1)*idim +3)=MOTS(3)   
             LISDUA(1+(km1-1)*idim +3)= NOMDU(ICOR(3))
             NOELEP(1+(km1-1)*idim +3)= KM1+1 
             NOELED(1+(km1-1)*idim +3)= KM1+1 
           endif
 347     continue
       endif
       
       SEGDES DESCR
       IRIGEL(3,INRI+I-1) = DESCR
       IRIGEL(5,INRI+I-1) = NIFOUR 
C       IRIGEL(6,INRI+I-1) = 0
328   CONTINUE       
c---- fin de boucle sur les xmatri ----
      SEGSUP MTRAV 

      IF (LISOUS(/1).NE.0) SEGDES IPT1    
395   CONTINUE
      SEGDES ICOUNT
C____ FIN DE BOUCLE SUR LES SOUS ZONES __________________________                                                
 
400   CONTINUE 


C**** MENAGE ET ECRITURE DE LA RIGIDTE RESULTATS **********************      
      SEGSUP ICOSOM
      SEGSUP IEXCLU 

      SEGDES MRIGID    
      SEGDES MELEME , IPT2 
      SEGDES MLMOTS
      CALL ECROBJ('RIGIDITE',MRIGID) 
 
      RETURN
      END