C FRIG3C    SOURCE    FANDEUR   11/10/07    21:15:23     7159           

      SUBROUTINE FRIG3C (maifro,IPRIGI,IPCHJE,IPRIG2)

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

*  ce sous-programme calcule la raideur de frottement en 3D
*  il a besoin pour cela du maillage de frottement et de la raideur
*  de contact (ou la raideur totale si c'est plus simple)

-INC CCOPTIO
-INC CCREEL

-INC SMCHPOI
-INC SMCOORD
-INC SMELEME
-INC SMMODEL
-INC SMRIGID

*  icpr lx du contact
      segment icpr(xcoor(/1)/idimp1)
      segment mwrk
       real*8 xdir(nbno),ydir(nbno),zdir(nbno)
       real*8 alpha(nbno),beta(nbno),gamma(nbno)
       integer in1(nbno),in2(nbno),in3(nbno),in4(nbno)
      endsegment
      segment mwrkf
       real*8 delta(nbno)
       integer in5(nbno),in6(nbno)
      endsegment

      character*4 modepl(3),moforc(3)
      data modepl /'UX  ','UY  ','UZ  '/
      data moforc /'FX  ','FY  ','FZ  '/
*
* test sur les options de calcul (faits dans frig)
C*      IF (IDIM .NE. 3) THEN
C*        CALL ERREUR(5)
C*        RETURN
C*      ENDIF
*
* quelques initialisations
*
      IPRIGI2 = 0
*
      idimp1 = IDIM + 1
*
      meleme = maifro
      segact,meleme
      if (lisous(/1).ne.0) then
        call erreur(16)
        GOTO 9900
      endif
      if (itypel.ne.22) call erreur(16)
      if (num(/1).ne.7 .and. num(/1).ne.9) call erreur(16)
      if (ierr.ne.0) GOTO 9900

* lecture de la raideur de contact et extraction de tous les LX
      mwrk = 0
      mwrkf = 0
      segini icpr

      nbno=0
      nbel=0
      nbel3=0
      nbel4=0
      nbel5=0
      nbel7=0
      xcr=xpetit*1.e+10

      mrigid = IPRIGI
      SEGACT,mrigid

      DO 10 irig = 1, irigel(/2)
        if (irigel(6,irig).eq.0) goto 10
        ipt2=irigel(1,irig)
        segact ipt2
        IF (ipt2.itypel.ne.22) goto 10
        ityp2  = ipt2.num(/1)
        nbel22 = ipt2.num(/2)
        nbel   = nbel + nbel22
c* Formulation forte/standard
        if (ityp2.EQ.3) THEN
          nbel3 = nbel3 + nbel22
        ELSEIF (ityp2.EQ.4) THEN
          nbel4 = nbel4 + nbel22
        ELSEIF (ityp2.EQ.5) THEN
          nbel5 = nbel5 + nbel22
c* Formulation faible
        ELSEIF (ityp2.EQ.7) THEN
          nbel7 = nbel7 + nbel22
        ELSE
          CALL ERREUR(16)
          goto 10
        ENDIF
        i1 = (ipt2.num(2,1)-1)*idimp1
        x1 = xcoor(i1+1)
        y1 = xcoor(i1+2)
        z1 = xcoor(i1+3)
        do 101 iel = 1, nbel22
          nbno=nbno+1
          icpr(ipt2.num(1,iel)) = nbno
          i2 = (ipt2.num(2,iel)-1)*idimp1
          xcr = MAX(xcr,ABS(x1-xcoor(i2+1)),
     &                  ABS(y1-xcoor(i2+2)),
     &                  ABS(z1-xcoor(i2+3)))
 101    continue
 10   CONTINUE

      xcr=xcr*1.E-12
*     write (ioimp,*) ' xcr ',xcr

      segini,mwrk
      if (nbel7.ne.0) segini,mwrkf

      DO 11 irig=1,irigel(/2)
        if (irigel(6,irig).eq.0) goto 11
        ipt2=irigel(1,irig)
        if (ipt2.itypel.ne.22) goto 110
        ityp2  = ipt2.num(/1)
        nbel22 = ipt2.num(/2)
        xmatri=irigel(4,irig)
        segact xmatri
c* Relation type 1 : noeud-noeud (forte/standard)
        if (ityp2.eq.3) then
          do 113 iel = 1, nbel22
            ip = icpr(ipt2.num(1,iel))
c*            if (ip.eq.0) goto 113
            in1(ip)=ipt2.num(2,iel)
            in2(ip)=ipt2.num(3,iel)
            xdir(ip)= re(1,2,iel)
            ydir(ip)= re(1,3,iel)
            zdir(ip)= re(1,4,iel)
c*            xdir(ip)=-re(1,5,iel)
c*            ydir(ip)=-re(1,6,iel)
c*            zdir(ip)=-re(1,7,iel)
c*            alpha(ip) = 1.
*        write(ioimp,*) ' ip 1 ',ip,xdir(ip),ydir(ip),zdir(ip)
 113      continue
c* Relation type 2 : noeud-segment (forte/standard)
        elseif (ityp2.eq.4) then
          do 114 iel = 1, nbel22
            ip = icpr(ipt2.num(1,iel))
c*            if (ip.eq.0) goto 114
            in1(ip)=ipt2.num(2,iel)
            in2(ip)=ipt2.num(3,iel)
            in3(ip)=ipt2.num(4,iel)
            xdir(ip)=-re(1, 8,iel)
            ydir(ip)=-re(1, 9,iel)
            zdir(ip)=-re(1,10,iel)
            if (abs(xdir(ip)).gt.0.1) then
              alpha(ip)= re(1,2,iel)/xdir(ip)
              beta(ip) = re(1,5,iel)/xdir(ip)
            elseif (abs(ydir(ip)).gt.0.1) then
              alpha(ip)= re(1,3,iel)/ydir(ip)
              beta(ip) = re(1,6,iel)/ydir(ip)
            else
              alpha(ip)= re(1,4,iel)/zdir(ip)
              beta(ip) = re(1,7,iel)/zdir(ip)
            endif
c*            beta(ip) = 1. - alpha(ip)
*        write(ioimp,*) ' ip 2 ',ip,xdir(ip),ydir(ip),zdir(ip),
*    >     alpha(ip),beta(ip)
 114      continue
c* Relation type 3 : noeud-triangle (forte/standard)
        elseif (ityp2.eq.5) then
          do 115 iel = 1, nbel22
            ip = icpr(ipt2.num(1,iel))
c*            if (ip.eq.0) goto 115
            in1(ip)=ipt2.num(2,iel)
            in2(ip)=ipt2.num(3,iel)
            in3(ip)=ipt2.num(4,iel)
            in4(ip)=ipt2.num(5,iel)
            xdir(ip)=-re(1,11,iel)
            ydir(ip)=-re(1,12,iel)
            zdir(ip)=-re(1,13,iel)
            if (abs(xdir(ip)).gt.0.1) then
              alpha(ip)= re(1,2,iel)/xdir(ip)
              beta(ip) = re(1,5,iel)/xdir(ip)
              gamma(ip)= re(1,8,iel)/xdir(ip)
            elseif (abs(ydir(ip)).gt.0.1) then
              alpha(ip)= re(1,3,iel)/ydir(ip)
              beta(ip) = re(1,6,iel)/ydir(ip)
              gamma(ip)= re(1,9,iel)/ydir(ip)
            else
              alpha(ip)= re(1, 4,iel)/zdir(ip)
              beta(ip) = re(1, 7,iel)/zdir(ip)
              gamma(ip)= re(1,10,iel)/zdir(ip)
            endif
c*            gamma(ip)= 1. - beta(ip) - gamma(ip)
*        write (ioimp,*) ' ip 3 ',ip,xdir(ip),ydir(ip),zdir(ip),
*    >     alpha(ip),beta(ip),gamma(ip)
 115      continue
c* Relation type 0 : triangle-triangle (faible)
        elseif (ityp2.eq.7) then
          do 117 iel = 1, nbel22
            ip = icpr(ipt2.num(1,iel))
c*            if (ip.eq.0) goto 117
            in1(ip)=ipt2.num(2,iel)
            in2(ip)=ipt2.num(3,iel)
            in3(ip)=ipt2.num(4,iel)
            in4(ip)=ipt2.num(5,iel)
            in5(ip)=ipt2.num(6,iel)
            in6(ip)=ipt2.num(7,iel)
            xdir(ip)=re(1, 2,iel)+re(1, 5,iel)+re(1, 8,iel)
            ydir(ip)=re(1, 3,iel)+re(1, 6,iel)+re(1, 9,iel)
            zdir(ip)=re(1, 4,iel)+re(1, 7,iel)+re(1,10,iel)
            if (abs(xdir(ip)).gt.0.1) then
              alpha(ip)= re(1, 2,iel)/xdir(ip)
              beta(ip) = re(1, 5,iel)/xdir(ip)
              gamma(ip)=-re(1,11,iel)/xdir(ip)
              delta(ip)=-re(1,14,iel)/xdir(ip)
            elseif (abs(ydir(ip)).gt.0.1) then
              alpha(ip)= re(1, 3,iel)/ydir(ip)
              beta(ip) = re(1, 6,iel)/ydir(ip)
              gamma(ip)=-re(1,12,iel)/ydir(ip)
              delta(ip)=-re(1,15,iel)/ydir(ip)
            else
              alpha(ip)= re(1, 4,iel)/zdir(ip)
              beta(ip) = re(1, 7,iel)/zdir(ip)
              gamma(ip)=-re(1,13,iel)/zdir(ip)
              delta(ip)=-re(1,16,iel)/zdir(ip)
            endif
*        write (ioimp,*) ' ip 0 ',ip,xdir(ip),ydir(ip),zdir(ip),
*    >     alpha(ip),beta(ip),gamma(ip),delta(ip)
 117      continue
        endif
        segdes xmatri
 110    continue
        segdes ipt2
  11  CONTINUE
      segdes mrigid
      IF (ierr.NE.0) GOTO 9910

* lecture du champ de jeux (pas de frottement si jeu non nul)
      mchpoi = IPCHJE
      segact mchpoi
      do 15 isoupo=1,ipchp(/1)
        msoupo=ipchp(isoupo)
        segact msoupo
        mpoval=ipoval
        segact mpoval
        if (vpocha(/2).ne.1) call erreur(16)
        ipt8=igeoc
        segact ipt8
        if (ipt8.num(/1).ne.vpocha(/2)) call erreur(16)
        do 16 i=1,vpocha(/1)
          if ((vpocha(i,1)).gt.xcr) then
            icpr(ipt8.num(1,i))=-abs(icpr(ipt8.num(1,i)))
*          write(ioimp,*) ' frig3c frottement diminué ',ipt8.num(1,i),
*     &                   vpocha(i,1)
          else
*          write(ioimp,*) ' frig3c frottement gardé   ',ipt8.num(1,i),
*     &                   vpocha(i,1)
          endif
  16    continue
        segdes ipt8,mpoval,msoupo
  15  continue
      segdes mchpoi
      IF (ierr.NE.0) GOTO 9910
*
*  creation des raideurs de frottement la ou il existe une raideur de
*  contact
*
      nrigel=0
      if (nbel3.ne.0) nrigel=nrigel+1
      if (nbel4.ne.0) nrigel=nrigel+1
      if (nbel5.ne.0) nrigel=nrigel+1
      if (nbel7.ne.0) nrigel=nrigel+1
*      write(ioimp,*) ' nrigel ',nrigel,nbel3,nbel4,nbel5,nbel7
*
      segini mrigid
*
      ichole=0
      imgeo1=0
      imgeo2=0
      isupeq=0
      iforig=ifomod
      mtymat='RIGIDITE'
* tout ce qu'on cree est du frottement et est symetrique
      do irig = 1, nrigel
        irigel(6,irig) = 2
        irigel(7,irig) = 0
        coerig(irig) = 1.D0
      enddo

      irig=0
*
*  Relation type 3 : noeud-triangle (forte/standard)
*  -----------------
      if (nbel5.ne.0) then
*
        irig=irig+1
        nelri5=nbel5*2
*
*     creation du meleme associe a la relation
*
        nbelem=nelri5
        nbnn  =5
        nbsous=0
        nbref =0
        segini ipt2
        ipt2.itypel=22
        irigel(1,irig)=ipt2
*
*     creation du descriptif commun a toutes les raideurs
*
        nligrp=13
        nligrd=nligrp
        segini descr
        lisinc(1)='LX  '
        lisdua(1)='FLX '
        noelep(1)=1
        noeled(1)=1
        do 501 i=2,nligrp,3
          lisinc(i  )=modepl(1)
          lisinc(i+1)=modepl(2)
          lisinc(i+2)=modepl(3)
          lisdua(i  )=moforc(1)
          lisdua(i+1)=moforc(2)
          lisdua(i+2)=moforc(3)
          noelep(i  )=(i+4)/3
          noelep(i+1)=noelep(i)
          noelep(i+2)=noelep(i)
          noeled(i  )=noelep(i)
          noeled(i+1)=noelep(i)
          noeled(i+2)=noelep(i)
 501    continue
        segdes,descr
        irigel(3,irig)=descr
*
*     creation du segment xmatri
*
        nelrig=nelri5
        segini xmatri
        irigel(4,irig)=xmatri
*
*  on boucle sur le maillage de frottement en creant les raideurs uniquement
*  la ou il y a une raideur de contact
*
        nelri5=0
*
        do 500 iel = 1, num(/2)
*
          ip=icpr(num(1,iel))
          if (ip.eq.0) goto 500
*  coefficient multiplicateur selon signe de icpr ("frottement diminue")
          if (ip.lt.0) then
            xcof=0.5d0
          else
            xcof=1.d0
          endif
          ip=ABS(ip)
*
*   on range le maillage sous-jacent dans le meleme
*
          ip4 =in4(ip)
c* test ci-dessous ?
          if (ip4.eq.0) goto 500
          ip3 =in3(ip)
          ip2 =in2(ip)
          ip1 =in1(ip)
          nelri5=nelri5+1
          ipt2.num(1,nelri5)=num(6,iel)
          ipt2.num(2,nelri5)=ip1
          ipt2.num(3,nelri5)=ip2
          ipt2.num(4,nelri5)=ip3
          ipt2.num(5,nelri5)=ip4
          melri5=nelri5
          nelri5=nelri5+1
          ipt2.num(1,nelri5)=num(7,iel)
          ipt2.num(2,nelri5)=ip1
          ipt2.num(3,nelri5)=ip2
          ipt2.num(4,nelri5)=ip3
          ipt2.num(5,nelri5)=ip4
*
*   on remplit les xmatri
*
*   directions orthogonales à xdir ydir zdir
*
          xr=-2.*zdir(ip)
          yr=-3.*zdir(ip)
          zr= 2.*xdir(ip)+3.*ydir(ip)
          if (ABS(xr).le.0.1.AND.ABS(yr).le.0.1.AND.ABS(zr).LE.0.1) then
            xr=-2.*ydir(ip)
            yr= 2.*xdir(ip)+3.*zdir(ip)
            zr=-3.*ydir(ip)
          endif
          sr=SQRT(xr*xr+yr*yr+zr*zr)
          xr=xr/sr
          yr=yr/sr
          zr=zr/sr
          xs=ydir(ip)*zr-zdir(ip)*yr
          ys=zdir(ip)*xr-xdir(ip)*zr
          zs=xdir(ip)*yr-ydir(ip)*xr
          ss=SQRT(xs*xs+ys*ys+zs*zs)
          xs=xs/ss
          ys=ys/ss
          zs=zs/ss
*
*  coefficient multiplicateur selon signe de icpr ("frottement diminue")
*  corrections sur les directions (xyz)r et (xyz)s
          alp_d = alpha(ip)*xcof
          bet_d =  beta(ip)*xcof
          gam_d = gamma(ip)*xcof
*
*  lambda
          re(1,1,melri5) = 0.
          re(1,1,nelri5) = 0.
*  ip1
          re(2,1,melri5) = alp_d*xr
          re(3,1,melri5) = alp_d*yr
          re(4,1,melri5) = alp_d*zr
          re(2,1,nelri5) = alp_d*xs
          re(3,1,nelri5) = alp_d*ys
          re(4,1,nelri5) = alp_d*zs
*  ip2
          re(5,1,melri5) = bet_d*xr
          re(6,1,melri5) = bet_d*yr
          re(7,1,melri5) = bet_d*zr
          re(5,1,nelri5) = bet_d*xs
          re(6,1,nelri5) = bet_d*ys
          re(7,1,nelri5) = bet_d*zs
*  ip3
          re( 8,1,melri5) = gam_d*xr
          re( 9,1,melri5) = gam_d*yr
          re(10,1,melri5) = gam_d*zr
          re( 8,1,nelri5) = gam_d*xs
          re( 9,1,nelri5) = gam_d*ys
          re(10,1,nelri5) = gam_d*zs
*  jp = ip4
          re(11,1,melri5) = -xr*xcof
          re(12,1,melri5) = -yr*xcof
          re(13,1,melri5) = -zr*xcof
          re(11,1,nelri5) = -xs*xcof
          re(12,1,nelri5) = -ys*xcof
          re(13,1,nelri5) = -zs*xcof
*  on transpose
          do 510 ic=2,nligrp
            re(1,ic,melri5)=re(ic,1,melri5)
            re(1,ic,nelri5)=re(ic,1,nelri5)
 510      continue
*
 500    continue
*
        if (nelrig.ne.nelri5) then
          nelrig=nelri5
          segadj,xmatri
          nbelem=nelri5
          nbnn  =5
          nbsous=0
          nbref =0
          segadj,ipt2
        endif
*
        segdes,xmatri,ipt2
      endif
*
*  Relation type 2 : noeud-segment (forte/standard)
*  -----------------
      if (nbel4.ne.0) then
*
        irig=irig+1
        nelri4 = nbel4*2

*     creation du meleme associe a la relation
*
        nbelem=nelri4
        nbnn  =4
        nbsous=0
        nbref =0
        segini ipt2
        ipt2.itypel=22
        irigel(1,irig)=ipt2
*
*     creation du descriptif commun a toutes les raideurs
*
        nligrp=10
        nligrd=nligrp
        segini,descr
        lisinc(1)='LX  '
        lisdua(1)='FLX '
        noelep(1)=1
        noeled(1)=1
        do 401 i=2,nligrp,3
          lisinc(i  )=modepl(1)
          lisinc(i+1)=modepl(2)
          lisinc(i+2)=modepl(3)
          lisdua(i  )=moforc(1)
          lisdua(i+1)=moforc(2)
          lisdua(i+2)=moforc(3)
          noelep(i  )=(i+4)/3
          noelep(i+1)=noelep(i)
          noelep(i+2)=noelep(i)
          noeled(i  )=noelep(i)
          noeled(i+1)=noelep(i)
          noeled(i+2)=noelep(i)
 401    continue
        segdes,descr
        irigel(3,irig)=descr
*
*     creation du segment xmatri
*
        nelrig = nelri4
        segini xmatri
        irigel(4,irig)=xmatri
*
*  on boucle sur le maillage de frottement en creant les raideurs uniquement
*  la ou il y a une raideur de contact
*
        nelri4=0
*
        do 400 iel=1,num(/2)
          ip=icpr(num(1,iel))
          if (ip.eq.0) goto 400
*  coefficient multiplicateur selon signe de ip ("frottement diminue")
          if (ip.lt.0) then
            xcof=0.5d0
          else
            xcof=1.d0
          endif
          ip=ABS(ip)
*
*   on range le maillage sous-jacent dans le meleme
*
c* Voir les tests ci-dessous ?
          if (in4(ip).ne.0) goto 400
          ip3 =in3(ip)
          if (ip3.eq.0) goto 400
          ip1 =in1(ip)
          ip2 =in2(ip)
          nelri4=nelri4+1
          ipt2.num(1,nelri4)=num(6,iel)
          ipt2.num(2,nelri4)=ip1
          ipt2.num(3,nelri4)=ip2
          ipt2.num(4,nelri4)=ip3
          melri4=nelri4
          nelri4=nelri4+1
          ipt2.num(1,nelri4)=num(7,iel)
          ipt2.num(2,nelri4)=ip1
          ipt2.num(3,nelri4)=ip2
          ipt2.num(4,nelri4)=ip3
*
*   on remplit les xmatri
*
*   directions orthogonales à xdir ydir zdir
*
          xr=-2.*zdir(ip)
          yr=-3.*zdir(ip)
          zr=2.*xdir(ip)+3.*ydir(ip)
          if (abs(xr).le.0.1.and.abs(yr).le.0.1.and.abs(zr).le.0.1) then
            xr=-2.*ydir(ip)
            yr=2.*xdir(ip)+3.*zdir(ip)
            zr=-3.*ydir(ip)
          endif
          sr=sqrt(xr*xr+yr*yr+zr*zr)
          xr=xr/sr
          yr=yr/sr
          zr=zr/sr
          xs=ydir(ip)*zr-zdir(ip)*yr
          ys=zdir(ip)*xr-xdir(ip)*zr
          zs=xdir(ip)*yr-ydir(ip)*xr
          ss=sqrt(xs*xs+ys*ys+zs*zs)
          xs=xs/ss
          ys=ys/ss
          zs=zs/ss
*
*  coefficient multiplicateur selon signe de icpr ("frottement diminue")
*  corrections sur les directions (xyz)r et (xyz)s
          alp_d = alpha(ip)*xcof
          bet_d =  beta(ip)*xcof
*  lambda
          re(1,1,melri4) = 0.
          re(1,1,nelri4) = 0.
*  ip1
          re(2,1,melri4) = alp_d*xr
          re(3,1,melri4) = alp_d*yr
          re(4,1,melri4) = alp_d*zr
          re(2,1,nelri4) = alp_d*xs
          re(3,1,nelri4) = alp_d*ys
          re(4,1,nelri4) = alp_d*zs
*  ip2
          re(5,1,melri4) = bet_d*xr
          re(6,1,melri4) = bet_d*yr
          re(7,1,melri4) = bet_d*zr
          re(5,1,nelri4) = bet_d*xs
          re(6,1,nelri4) = bet_d*ys
          re(7,1,nelri4) = bet_d*zs
*  jp = ip3
          re( 8,1,melri4) = -xr*xcof
          re( 9,1,melri4) = -yr*xcof
          re(10,1,melri4) = -zr*xcof
          re( 8,1,nelri4) = -xs*xcof
          re( 9,1,nelri4) = -ys*xcof
          re(10,1,nelri4) = -zs*xcof
*
*  on transpose
          do 410 ic=2,nligrp
            re(1,ic,melri4)=re(ic,1,melri4)
            re(1,ic,nelri4)=re(ic,1,nelri4)
 410      continue
 400    continue
*
        if (nelrig.ne.nelri4) then
          nelrig=nelri4
          segadj,xmatri
          nbelem=nelri4
          nbnn  =4
          nbsous=0
          nbref =0
          segadj,ipt2
        endif
*
        segdes xmatri,ipt2,descr
      endif
*
*  Relation type 1 : noeud-noeud (forte/standard)
*  -----------------
      if (nbel3.ne.0) then
*
        irig=irig+1
        nelri3 = nbel3*2
*
*     creation du meleme associe a la relation
*
        nbelem = nelri3
        nbnn   = 3
        nbsous = 0
        nbref  = 0
        segini ipt2
        ipt2.itypel=22
        irigel(1,irig)=ipt2
*
*     creation du descriptif commun a toutes les raideurs
*
        nligrp=7
        nligrd=nligrp
        segini,descr
        lisinc(1)='LX  '
        lisdua(1)='FLX '
        noelep(1)=1
        noeled(1)=1
        do 301 i=2,nligrp,3
          lisinc(i  )=modepl(1)
          lisinc(i+1)=modepl(2)
          lisinc(i+2)=modepl(3)
          lisdua(i  )=moforc(1)
          lisdua(i+1)=moforc(2)
          lisdua(i+2)=moforc(3)
          noelep(i  )=(i+4)/3
          noelep(i+1)=noelep(i)
          noelep(i+2)=noelep(i)
          noeled(i  )=noelep(i)
          noeled(i+1)=noelep(i)
          noeled(i+2)=noelep(i)
 301    continue
        segdes,descr
        irigel(3,irig)=descr
*
*     creation du segment xmatri
*
        nelrig = nelri3
        segini xmatri
        irigel(4,irig)=xmatri
*
*  on boucle sur le maillage de frottement en creant les raideurs uniquement
*  la ou il y a une raideur de contact
*
        nelri3=0
*
        do 300 iel=1,num(/2)
          ip=icpr(num(1,iel))
          if (ip.eq.0) goto 300
*  coefficient multiplicateur selon signe de ip ("frottement diminue")
          if (ip.lt.0) then
            xcof=0.5d0
          else
            xcof=1.d0
          endif
          ip=ABS(ip)
*
*   on range le maillage sous-jacent dans le meleme
*
c* Voir les tests ci-dessous ?
          if (in3(ip).ne.0) goto 300
          if (in4(ip).ne.0) goto 300
          ip1 =in1(ip)
          ip2 =in2(ip)
          nelri3=nelri3+1
          ipt2.num(1,nelri3)=num(6,iel)
          ipt2.num(2,nelri3)=ip1
          ipt2.num(3,nelri3)=ip2
          melri3=nelri3
          nelri3=nelri3+1
          ipt2.num(1,nelri3)=num(7,iel)
          ipt2.num(2,nelri3)=ip1
          ipt2.num(3,nelri3)=ip2
*
*   on remplit les xmatri
*
*   directions orthogonales à xdir ydir zdir
*
          xr=-2.*zdir(ip)
          yr=-3.*zdir(ip)
          zr=2.*xdir(ip)+3.*ydir(ip)
          if (abs(xr).le.0.1.and.abs(yr).le.0.1.and.abs(zr).le.0.1) then
            xr=-2.*ydir(ip)
            yr=2.*xdir(ip)+3.*zdir(ip)
            zr=-3.*ydir(ip)
          endif
          sr=sqrt(xr*xr+yr*yr+zr*zr)
          xr=xr/sr
          yr=yr/sr
          zr=zr/sr
          xs=ydir(ip)*zr-zdir(ip)*yr
          ys=zdir(ip)*xr-xdir(ip)*zr
          zs=xdir(ip)*yr-ydir(ip)*xr
          ss=sqrt(xs*xs+ys*ys+zs*zs)
          xs=xs/ss
          ys=ys/ss
          zs=zs/ss
*
*  coefficient multiplicateur selon signe de icpr ("frottement diminue")
*  corrections sur les directions (xyz)r et (xyz)s
          xr = xr*xcof
          yr = yr*xcof
          zr = zr*xcof
          xs = xs*xcof
          ys = ys*xcof
          zs = zs*xcof
*  lambda
          re(1,1,melri3)= 0.
          re(1,1,nelri3)= 0.
*  ip1
          re(2,1,melri3)= xr
          re(3,1,melri3)= yr
          re(4,1,melri3)= zr
          re(2,1,nelri3)= xs
          re(3,1,nelri3)= ys
          re(4,1,nelri3)= zs
*  jp = ip2
          re(5,1,melri3)= -xr
          re(6,1,melri3)= -yr
          re(7,1,melri3)= -zr
          re(5,1,nelri3)= -xs
          re(6,1,nelri3)= -ys
          re(7,1,nelri3)= -zs
*  on transpose
          do 310 ic=2,nligrp
            re(1,ic,melri3)=re(ic,1,melri3)
            re(1,ic,nelri3)=re(ic,1,nelri3)
 310      continue
*
 300    continue
*
        if (nelrig.ne.nelri3) then
          nelrig=nelri3
          segadj,xmatri
          nbelem=nelri3
          nbnn  =3
          nbsous=0
          nbref =0
          segadj,ipt2
        endif
*
        segdes xmatri,ipt2
      endif
*
*  Relation type 0 : triangle-triangle (faible)
*  -----------------
      if (nbel7.ne.0) then
*
        irig=irig+1
        nelri7=nbel7*2
*
*     creation du meleme associe a la relation
*
        nbelem=nelri7
        nbnn  =7
        nbsous=0
        nbref =0
        segini ipt2
        ipt2.itypel=22
        irigel(1,irig)=ipt2
*
*     creation du descriptif commun a toutes les raideurs
*        nligrp=idim*(nbnn-1)+1
*
        nligrp=19
        nligrd=nligrp
        segini,descr
        lisinc(1)='LX  '
        lisdua(1)='FLX '
        noelep(1)=1
        noeled(1)=1
        do 701 i=2,nligrp,3
          lisinc(i  )=modepl(1)
          lisinc(i+1)=modepl(2)
          lisinc(i+2)=modepl(3)
          lisdua(i  )=moforc(1)
          lisdua(i+1)=moforc(2)
          lisdua(i+2)=moforc(3)
          noelep(i  )=(i+4)/3
          noelep(i+1)=noelep(i)
          noelep(i+2)=noelep(i)
          noeled(i  )=noelep(i)
          noeled(i+1)=noelep(i)
          noeled(i+2)=noelep(i)
 701    continue
        segdes,descr
        irigel(3,irig)=descr
*
*     creation du segment xmatri
*
        nelrig=nelri7
        segini xmatri
        irigel(4,irig)=xmatri
*
*  on boucle sur le maillage de frottement en creant les raideurs uniquement
*  la ou il y a une raideur de contact
*
        nelri7=0
        do 700 iel = 1, num(/2)
          ip=icpr(num(1,iel))
          if (ip.eq.0) goto 700
*  coefficient multiplicateur selon signe de icpr ("frottement diminue")
          if (ip.lt.0) then
            xcof=0.5d0
          else
            xcof=1.d0
          endif
          ip=ABS(ip)
*
*   on range le maillage sous-jacent dans le meleme
*
c* Voir les tests ci-dessous ?
          ip4 = in4(ip)
          if (ip4.eq.0) goto 700
          ip5 = in5(ip)
          if (ip5.eq.0) goto 700
          ip6 = in6(ip)
          if (ip6.eq.0) goto 700
          ip1 = in1(ip)
          ip2 = in2(ip)
          ip3 = in3(ip)
          nelri7=nelri7+1
          ipt2.num(1,nelri7)=num(8,iel)
          ipt2.num(2,nelri7)=ip1
          ipt2.num(3,nelri7)=ip2
          ipt2.num(4,nelri7)=ip3
          ipt2.num(5,nelri7)=ip4
          ipt2.num(6,nelri7)=ip5
          ipt2.num(7,nelri7)=ip6
          melri7=nelri7
          nelri7=nelri7+1
          ipt2.num(1,nelri7)=num(9,iel)
          ipt2.num(2,nelri7)=ip1
          ipt2.num(3,nelri7)=ip2
          ipt2.num(4,nelri7)=ip3
          ipt2.num(5,nelri7)=ip4
          ipt2.num(6,nelri7)=ip5
          ipt2.num(7,nelri7)=ip6
*
*   on remplit les xmatri
*
*   directions orthogonales à xdir ydir zdir
*
          xr=-2.*zdir(ip)
          yr=-3.*zdir(ip)
          zr= 2.*xdir(ip)+3.*ydir(ip)
          if (ABS(xr).le.0.1.AND.ABS(yr).le.0.1.AND.ABS(zr).LE.0.1) then
            xr=-2.*ydir(ip)
            yr= 2.*xdir(ip)+3.*zdir(ip)
            zr=-3.*ydir(ip)
          endif
          sr=SQRT(xr*xr+yr*yr+zr*zr)
          xr=xr/sr
          yr=yr/sr
          zr=zr/sr
          xs=ydir(ip)*zr-zdir(ip)*yr
          ys=zdir(ip)*xr-xdir(ip)*zr
          zs=xdir(ip)*yr-ydir(ip)*xr
          ss=SQRT(xs*xs+ys*ys+zs*zs)
          xs=xs/ss
          ys=ys/ss
          zs=zs/ss
*
*  coefficient multiplicateur selon signe de icpr ("frottement diminue")
*  corrections sur les directions (xyz)r et (xyz)s
          alp_d = alpha(ip)*xcof
          bet_d =  beta(ip)*xcof
          uab_d =  xcof - alp_d - bet_d
          gam_d = -gamma(ip)*xcof
          del_d = -delta(ip)*xcof
          ugd_d = -(xcof + gam_d + del_d)
*
*  lambda
          re( 1,1,melri7) = 0.
          re( 1,1,nelri7) = 0.
*  ip1
          re( 2,1,melri7) = alp_d*xr
          re( 3,1,melri7) = alp_d*yr
          re( 4,1,melri7) = alp_d*zr
          re( 2,1,nelri7) = alp_d*xs
          re( 3,1,nelri7) = alp_d*ys
          re( 4,1,nelri7) = alp_d*zs
*  ip2
          re( 5,1,melri7) = bet_d*xr
          re( 6,1,melri7) = bet_d*yr
          re( 7,1,melri7) = bet_d*zr
          re( 5,1,nelri7) = bet_d*xs
          re( 6,1,nelri7) = bet_d*ys
          re( 7,1,nelri7) = bet_d*zs
*  ip3
          re( 8,1,melri7) = uab_d*xr
          re( 9,1,melri7) = uab_d*yr
          re(10,1,melri7) = uab_d*zr
          re( 8,1,nelri7) = uab_d*xs
          re( 9,1,nelri7) = uab_d*ys
          re(10,1,nelri7) = uab_d*zs
*  ip4 = jp1
          re(11,1,melri7) = gam_d*xr
          re(12,1,melri7) = gam_d*yr
          re(13,1,melri7) = gam_d*zr
          re(11,1,nelri7) = gam_d*xs
          re(12,1,nelri7) = gam_d*ys
          re(13,1,nelri7) = gam_d*zs
*  ip5 = jp2
          re(14,1,melri7) = del_d*xr
          re(15,1,melri7) = del_d*yr
          re(16,1,melri7) = del_d*zr
          re(14,1,nelri7) = del_d*xs
          re(15,1,nelri7) = del_d*ys
          re(16,1,nelri7) = del_d*zs
*  ip6 = jp3
          re(17,1,melri7) = ugd_d*xr
          re(18,1,melri7) = ugd_d*yr
          re(19,1,melri7) = ugd_d*zr
          re(17,1,nelri7) = ugd_d*xs
          re(18,1,nelri7) = ugd_d*ys
          re(19,1,nelri7) = ugd_d*zs
*  on transpose
          do 710 ic=2,nligrp
            re(1,ic,melri7)=re(ic,1,melri7)
            re(1,ic,nelri7)=re(ic,1,nelri7)
 710      continue
*
 700    continue
*
        if (nelrig.ne.nelri7) then
          nelrig=nelri7
          segadj,xmatri
          nbelem=nelri7
          nbnn  =7
          nbsous=0
          nbref =0
          segadj,ipt2
        endif
        segdes,xmatri,ipt2
      endif

C*      SEGDES,mrigid
      IPRIG2 = mrigid

 9910 CONTINUE
      SEGSUP,icpr,mwrk
      if (mwrkf.ne.0) segsup,mwrkf
 9900 CONTINUE
      SEGDES,meleme

      RETURN
      END