C FRIG2C    SOURCE    FANDEUR   11/10/07    21:15:21     7159           

      SUBROUTINE FRIG2C (maifro,IPRIGI,IPCHJE,IPRIG2)

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

*  Ce sous-programme calcule la raideur de frottement en 2D.
*  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 SMRIGID

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

      character*4 modepl(4),moforc(4)
      data modepl /'UX  ','UY  ','UR  ','UZ  '/
      data moforc /'FX  ','FY  ','FR  ','FZ  '/
*
* test sur les options de calcul (faits dans frig)
C*      IF (IDIM .NE. 2) THEN
C*        CALL ERREUR(5)
C*        RETURN
C*      ENDIF
C*      IF (IFOMOD .NE. -1 .AND. IFOMOD .NE. 0) THEN
C*        CALL ERREUR(710)
C*        RETURN
C*      ENDIF
*
* quelques initialisations
*
      IPRIG2 = 0
      idimp1 = IDIM + 1
*     indice servant a affecter le nom des composantes
      idepl = (ifomod + 1)*2
*
      meleme = maifro
      SEGACT,meleme
      if (lisous(/1).ne.0) then
        call erreur(16)
      ELSE IF (itypel.NE.22) then
        CALL ERREUR(16)
      ELSE IF (num(/1).NE.5 .AND. num(/1).NE.6) THEN
        CALL ERREUR(16)
      ENDIF
      IF (IERR.NE.0) GOTO 9900

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

      nbno=0
      nbel=0
      nbel3=0
      nbel4=0
      nbel5=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
c* Formulation faible
        ELSEIF (ityp2.EQ.5) THEN
          nbel5 = nbel5 + nbel22
        ELSE
          CALL ERREUR(16)
          goto 10
        ENDIF
        i1 = (ipt2.num(2,1)-1)*idimp1
        x1 = xcoor(i1+1)
        y1 = xcoor(i1+2)
        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)))
 101    continue
 10   CONTINUE

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

      segini mwrk

      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,4,iel)
            ydir(ip)=re(1,5,iel)
c*            alpha(ip)=+1.
c*            write(ioimp,*) ' ip            ',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,6,iel)
            ydir(ip)=re(1,7,iel)
            if (ABS(xdir(ip)).gt.0.1) then
              alpha(ip)=-re(1,2,iel)/xdir(ip)
            else
              alpha(ip)=-re(1,3,iel)/ydir(ip)
            endif
            beta(ip) = 1. - alpha(ip)
c*            write(ioimp,*) ' ip alpha      ',ip,alpha(ip)
 114      continue
c* Relation type 0 : segment-segment (faible)
        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,2,iel)+re(1,4,iel)
            ydir(ip)=re(1,3,iel)+re(1,5,iel)
* Normale de norme = 1
*           xl=sqrt(xdir(ip)**2+ydir(ip)**2)
*           if (ABS(xdir(ip)/xl).gt.0.1) then
            if (ABS(xdir(ip)).gt.0.1) then
              alpha(ip)= re(1,2,iel)/xdir(ip)
              beta(ip) =-re(1,6,iel)/xdir(ip)
            else
              alpha(ip)= re(1,3,iel)/ydir(ip)
              beta(ip) =-re(1,7,iel)/ydir(ip)
            endif
c*            write(ioimp,*) ' ip alpha beta ',ip,alpha(ip),beta(ip)
 115      continue
        endif
        segdes xmatri
 110    continue
        segdes ipt2
 11   CONTINUE
      segdes mrigid
      IF (ierr.NE.0) GOTO 9910

* Traitement du champ de jeux (demi-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,*) ' frig2c frottement diminué ',ipt8.num(1,i),
*    &     vpocha(i,1)
          else
*     write(ioimp,*) ' frig2c 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
*     write(ioimp,*) ' nrigel ',nrigel,nbel3,nbel4,nbel5
*
      segini mrigid
*
*  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
      ichole=0
      imgeo1=0
      imgeo2=0
      isupeq=0
      iforig=ifomod
      mtymat='RIGIDITE'
*
      irig=0
*
* Relation type 0 : segment-segment (faible)
*------------------
      if (nbel5.ne.0) then
        irig=irig+1
        nelri5=nbel5
*
*    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=9
        nligrd=nligrp
        segini descr
        lisinc(1)='LX  '
        lisdua(1)='FLX '
        noelep(1)=1
        noeled(1)=1
        do 501 i=2,9,2
          lisinc(i  )=modepl(1+idepl)
          lisinc(i+1)=modepl(2+idepl)
          lisdua(i  )=moforc(1+idepl)
          lisdua(i+1)=moforc(2+idepl)
          noelep(i  )=(i+2)/2
          noelep(i+1)=noelep(i)
          noeled(i  )=noelep(i)
          noeled(i+1)=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 la raideur 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)
          if (ip4.eq.0) goto 500
          ip1 =in1(ip)
          ip2 =in2(ip)
          ip3 =in3(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
*
          xr =  xcof*ydir(ip)
          yr = -xcof*xdir(ip)
*
*   on remplit le xmatri
*
*  lambda
          re(1,1,nelri5) = 0.
*  ip1
          re(2,1,nelri5) = alpha(ip)*xr
          re(3,1,nelri5) = alpha(ip)*yr
*  ip2
          re(4,1,nelri5) = (1.-alpha(ip))*xr
          re(5,1,nelri5) = (1.-alpha(ip))*yr
*  ip3
          re(6,1,nelri5) = -beta(ip)*xr
          re(7,1,nelri5) = -beta(ip)*yr
*  ip
          re(8,1,nelri5) = (beta(ip)-1.)*xr
          re(9,1,nelri5) = (beta(ip)-1.)*yr
*  on transpose
          do 510 ic=2,nligrd
            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
*
*     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=7
        nligrd=nligrp
        segini,descr
        lisinc(1)='LX  '
        lisdua(1)='FLX '
        noelep(1)=1
        noeled(1)=1
        do 401 i=2,nligrp,2
          lisinc(i)  =modepl(1+idepl)
          lisinc(i+1)=modepl(2+idepl)
          lisdua(i)  =moforc(1+idepl)
          lisdua(i+1)=moforc(2+idepl)
          noelep(i)  =(i+2)/2
          noelep(i+1)=noelep(i)
          noeled(i)  =noelep(i)
          noeled(i+1)=noelep(i)
 401    continue
        segdes,descr
        irigel(3,irig)=descr
*
*     creation du segment xmatri
*
        nelrig=nbel4
        segini xmatri
        irigel(4,irig)=xmatri
*
*  on boucle sur le maillage de frottement en creant la raideur 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 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* Tests a voir
          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(5,iel)
          ipt2.num(2,nelri4)=ip1
          ipt2.num(3,nelri4)=ip2
          ipt2.num(4,nelri4)=ip3
*
          xr =  xcof*ydir(ip)
          yr = -xcof*xdir(ip)
c* OF A voir !!!!!!!!!!!!!!!!!!!
c*          xr = -xcof*ydir(ip)
c*          yr =  xcof*xdir(ip)
*
*   on remplit le xmatri
*
*  lambda
          re(1,1,nelri4) = 0.
*  ip1
          re(2,1,nelri4) = alpha(ip)*xr
          re(3,1,nelri4) = alpha(ip)*yr
*  ip2
          re(4,1,nelri4) = beta(ip)*xr
          re(5,1,nelri4) = beta(ip)*yr
*  ip
          re(6,1,nelri4) = -xr
          re(7,1,nelri4) = -yr
*  on transpose
          do 410 ic=2,nligrp
            re(1,ic,nelri4)=re(ic,1,nelri4)
 410      continue
*
 400    continue

        if (nelrig.ne.nelri4) then
          nelrig=nelri4
          segadj,xmatri
          nbelem=nelri4
          nbsous=0
          nbref=0
          nbnn=4
          segadj,ipt2
        endif
        segdes xmatri,ipt2
      endif
*
*  Relation type 1 : noeud-noeud (forte/standard)
*
      if (nbel3.ne.0) then
*
        irig=irig+1
        nelri3=nbel3
*
*     creation du meleme associe a la relation
*
        nbsous=0
        nbref=0
        nbnn=3
        nbelem=nelri3
        segini ipt2
        ipt2.itypel=22
        irigel(1,irig)=ipt2
*
*     creation du descriptif commun a toutes les raideurs
*
        nligrp=5
        nligrd=nligrp
        segini descr
        lisinc(1)='LX  '
        lisdua(1)='FLX '
        noelep(1)=1
        noeled(1)=1
        do 301 i=2,nligrp,2
          lisinc(i)  =modepl(1+idepl)
          lisinc(i+1)=modepl(2+idepl)
          lisdua(i)  =moforc(1+idepl)
          lisdua(i+1)=moforc(2+idepl)
          noelep(i)  =(i+2)/2
          noelep(i+1)=noelep(i)
          noeled(i)  =noelep(i)
          noeled(i+1)=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 la raideur 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 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
*
          if (in3(ip).ne.0) goto 300
          ip1 =in1(ip)
          ip2 =in2(ip)
          nelri3=nelri3+1
          ipt2.num(1,nelri3)=num(5,iel)
          ipt2.num(2,nelri3)=ip1
          ipt2.num(3,nelri3)=ip2
*
          xr =  xcof*ydir(ip)
          yr = -xcof*xdir(ip)
*
*   on remplit le xmatri
*
*  lambda
          re(1,1,nelri3)= 0.
*  ip1
          re(2,1,nelri3)=  xr
          re(3,1,nelri3)=  yr
*  ip
          re(4,1,nelri3)= -xr
          re(5,1,nelri3)= -yr
*  on transpose
          do 310 ic=2,nligrd
            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

*     segdes,mrigid
      IPRIG2 = mrigid

 9910 CONTINUE
      SEGSUP,icpr,mwrk
 9900 CONTINUE
      SEGDES,meleme

      RETURN
      END