Télécharger impo32.eso

Retour à la liste

impo32
C IMPO32    SOURCE    CB215821  25/04/08    21:15:21     12227          
 
*  impo bloc en 3D
 
      SUBROUTINE IMPO32(ipt1,ipt6,ipt8,itcont,mchel1,mrigid,mchpoi)
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
-INC CCGEOME
 
-INC SMCOORD
-INC SMELEME
-INC SMRIGID
-INC SMCHPOI
-INC SMCHAML
 
* Contact en formulation forte :
* Segments utiles dans les cas dits pathologiques
      segment mfopa1
        integer lpoin(nfopa1,3)
        real*8  xpoin(nfopa1,3)
        real*8  zje1(nfopa1)
      endsegment
      segment mfopa2
        integer lsegm(nfopa2,4)
        real*8  xsegm(nfopa2,3)
        real*8  zje2(nfopa2)
      endsegment
* directions moyennes aux sommets
      segment mfopa3
        real*8 xms(nbpts),yms(nbpts),zms(nbpts)
        integer ims(nbpts)
      endsegment
* positions compressees
      segment mfopa4
        real*8 tco(nbel1),xw(nbel1)
        integer ipnum(nbel1),ico(nbel1),iw(nbel1)
      endsegment
      segment mfopa5
        integer ind(indt)
      endsegment
 
* Contact en formulation faible :
* Segment utile pour le calcul de l'intersection des 2 triangles
* projetes dans le plan de contact intermediaire
      segment mfaible
        real*8 cPT1C(3,3),cPT2C(3,3), cPT1A(4,2),cPT2A(4,2)
        real*8 cPIn0(6,2), cPIn(6,2), test(6)
        real*8 vT1A(3,2), vT2A(3,2)
        real*8 SuT1A, SuT2A, SuPIn
*DBG-F        real*8 xGIn,yGIn, b1T1,b2T1,b3T1, b1T2,b2T2,b3T2
      endsegment
 
      PARAMETER ( X1s3=0.333333333333333333333333333333333333333333D0 )
      PARAMETER ( sqr3=1.7320508075688772935D0)
 
      character*4 modepl(3),moforc(3)
 
      data modepl /'UX  ','UY  ','UZ  '/
      data moforc /'FX  ','FY  ','FZ  '/
*
      idimp1 = IDIM + 1
*
*  Lecture du maillage support des conditions de contact
*  il s(agit la du premier maillage de contact, de type tri3
*
**    write(6,*) ' ipt1 dans impo32 ',ipt1
**    write(6,*) ' ipt6 dans impo32 ',ipt6
**    write(6,*) ' ipt8 dans impo32 ',ipt8
*
*  Activation du maillage et petites verifications
*
      segact ipt8
      segact ipt6
      segact ipt1
**    write(6,*) ' ipt1 2  dans impo32 ',ipt1
      nbno6 = ipt6.num(/1)
      nbel6 = ipt6.num(/2)
      nbno8 = ipt8.num(/1)
      nbel8 = ipt8.num(/2)
      nbel1 = ipt1.num(/2)
**    write(6,*) 'nbel6 ',nbel6
**    write(6,*) 'nbel1 ',nbel1
      if (ipt1.lisous(/1).ne.0) call erreur(25)
      if (ipt6.itypel.ne.4) call erreur(16)
      if (nbno6.ne.3 ) call erreur(16)
      if (ierr.ne.0) goto 900
*
*   Increment sur le nombre d'elements de contact retenus et utilises
*   dans le chpoint (mpoval et igeoc) et rigidite pour agrandir les
*   longueurs des segments adequats
      incnel = 2000
*
*  Creation de la raideur des conditions de contact
*  Remplissage de l'entete commun
*
      nrigel = 1
      segini,mrigid
      ichole = 0
      imgeo1 = 0
      imgeo2 = 0
      isupeq = 0
      iforig = ifour
      coerig(1) = 1.
      mtymat='RIGIDITE'
*
*     MCHAML materiau associe au MMODEL
      MELVA1 = 0
      MELVA2 = 0
      IF (MCHEL1.NE.0) THEN
        SEGACT, MCHEL1
*  recherche rapide d'un maillage correspondant dans le mchaml
        DO 210 n2 = 1,mchel1.imache(/1)
*        write(ioimp,*) ' n2 imache(n2) ipt1 ',n2,mchel1.imache(n2),ipt1
          if (mchel1.imache(n2).eq.ipt1) goto 220
 210    continue
 
        goto 230
 220    continue
        MCHAM1 = MCHEL1.ICHAML(n2)
        SEGACT, MCHAM1
        NBCO = MCHAM1.NOMCHE(/2)
        DO JCMP = 1,NBCO
          IF (MCHAM1.NOMCHE(JCMP).EQ.'JEU') THEN
            MELVA1 = MCHAM1.IELVAL(JCMP)
            SEGACT, MELVA1
            NELJ = MELVA1.VELCHE(/2)
            NPTELJ = min(MELVA1.VELCHE(/1),4)
C
C           Utilise pour le zonage
            xjmax = -xgrand
            xjmin =  xgrand
            DO IJI2=1,NELJ
              DO IJI1=1,MELVA1.VELCHE(/1) 
                xjmax=max(xjmax,MELVA1.VELCHE(IJI1,IJI2))
                xjmin=min(xjmin,MELVA1.VELCHE(IJI1,IJI2))
              ENDDO
            ENDDO
C
          ENDIF
          IF (MCHAM1.NOMCHE(JCMP).EQ.'ADHE') THEN
            MELVA2 = MCHAM1.IELVAL(JCMP)
            SEGACT, MELVA2
            NELA = MELVA2.VELCHE(/2)
            NPTELA = min(MELVA2.VELCHE(/1),4)
          ENDIF
        ENDDO
      ENDIF
 230  continue
*
*  Creation du chpoint de depi
*
      nat=1
      nsoupo=1
      segini mchpoi
      mtypoi='DEPIMP'
      mochde='engendré par impose'
      ifopoi=ifour
      jattri(1)=2
*
      nc=2
      IF (MELVA2.NE.0) nc=3
      segini msoupo
      ipchp(1)=msoupo
      nocomp(1)='FLX '
      nocomp(2)='SCAL'
      IF (MELVA2.NE.0) THEN
        nocomp(3)='FADH'
      ENDIF
*
      nbnn  =1
      nbelem=incnel
      nbref =0
      nbsous=0
      segini ipt7
      ipt7.itypel=1
      igeoc=ipt7
*
      n=incnel
      segini mpoval
      ipoval=mpoval
*
*  Nombre de noeuds dans le chpoint (en totalite) : ipt7 et mpoval
      nelch = 0
 
*=======================================================================
* Formulation "forte" (standard) du contact :
*=======================================================================
* Element de contact 3D a 3 noeuds
**    write(6,*) ' itcont dans impo32 ',itcont
*  itcont 2: formulation faible
      if (itcont.eq.2) goto 500
 
*  Calcul des directions moyennes de la boite d'encadrement, et de la taille max
      segact mcoord
      segini mfopa3
 
      xmin=xgrand
      xmax=-xgrand
      ymin=xgrand
      ymax=-xgrand
      zmin=xgrand
      zmax=-xgrand
      tamax=0.d0
 
      do 820 iel6=1,nbel6
*
        ip1=ipt6.num(1,iel6)
        ip2=ipt6.num(2,iel6)
        ip3=ipt6.num(3,iel6)
        ipv1 = (ip1-1)*idimp1
        xp1 = xcoor(ipv1+1)
        yp1 = xcoor(ipv1+2)
        zp1 = xcoor(ipv1+3)
        ipv2 = (ip2-1)*idimp1
        xp2 = xcoor(ipv2+1)
        yp2 = xcoor(ipv2+2)
        zp2 = xcoor(ipv2+3)
        ipv3 = (ip3-1)*idimp1
        xp3 = xcoor(ipv3+1)
        yp3 = xcoor(ipv3+2)
        zp3 = xcoor(ipv3+3)
 
        xmin=min(xmin,xp1,xp2,xp3)
        xmax=max(xmax,xp1,xp2,xp3)
        ymin=min(ymin,yp1,yp2,yp3)
        ymax=max(ymax,yp1,yp2,yp3)
        zmin=min(zmin,zp1,zp2,zp3)
        zmax=max(zmax,zp1,zp2,zp3)
 
*
*    normale au plan (123)
*
        x12 = xp2 - xp1
        y12 = yp2 - yp1
        z12 = zp2 - zp1
        x23 = xp3 - xp2
        y23 = yp3 - yp2
        z23 = zp3 - zp2
        xn = y12*z23 - z12*y23
        yn = z12*x23 - x12*z23
        zn = x12*y23 - y12*x23
        sn = xn*xn + yn*yn + zn*zn
        sn = max(sqrt(abs(sn)),xpetit*1d10)
        tamax=max(sn,tamax)
 
        xn = xn/sn
        yn = yn/sn
        zn = zn/sn
        ims(ip1)=1
        xms(ip1)=xms(ip1)+xn
        yms(ip1)=yms(ip1)+yn
        zms(ip1)=zms(ip1)+zn
        ims(ip2)=1
        xms(ip2)=xms(ip2)+xn
        yms(ip2)=yms(ip2)+yn
        zms(ip2)=zms(ip2)+zn
        ims(ip3)=1
        xms(ip3)=xms(ip3)+xn
        yms(ip3)=yms(ip3)+yn
        zms(ip3)=zms(ip3)+zn
 820  continue
**    write(6,*) 'xmin xmax ymin ymax zmin zmax',
**   >            xmin,xmax,ymin,ymax,zmin,zmax   
      tamax = tamax * 5.
 
      do 822 iel8=1,nbel8
*
        ip1=ipt8.num(1,iel8)
        ip2=ipt8.num(2,iel8)
        ip3=ipt8.num(3,iel8)
        ipv1 = (ip1-1)*idimp1
        xp1 = xcoor(ipv1+1)
        yp1 = xcoor(ipv1+2)
        zp1 = xcoor(ipv1+3)
        ipv2 = (ip2-1)*idimp1
        xp2 = xcoor(ipv2+1)
        yp2 = xcoor(ipv2+2)
        zp2 = xcoor(ipv2+3)
        ipv3 = (ip3-1)*idimp1
        xp3 = xcoor(ipv3+1)
        yp3 = xcoor(ipv3+2)
        zp3 = xcoor(ipv3+3)
*
*    normale au plan (123)
*
        x12 = xp2 - xp1
        y12 = yp2 - yp1
        z12 = zp2 - zp1
        x23 = xp3 - xp2
        y23 = yp3 - yp2
        z23 = zp3 - zp2
        xn = y12*z23 - z12*y23
        yn = z12*x23 - x12*z23
        zn = x12*y23 - y12*x23
        sn = xn*xn + yn*yn + zn*zn
        sn = max(sqrt(abs(sn)),xpetit*1d10)
        xn = xn/sn
        yn = yn/sn
        zn = zn/sn
        ims(ip1)=1
        xms(ip1)=xms(ip1)+xn
        yms(ip1)=yms(ip1)+yn
        zms(ip1)=zms(ip1)+zn
        ims(ip2)=1
        xms(ip2)=xms(ip2)+xn
        yms(ip2)=yms(ip2)+yn
        zms(ip2)=zms(ip2)+zn
        ims(ip3)=1
        xms(ip3)=xms(ip3)+xn
        yms(ip3)=yms(ip3)+yn
        zms(ip3)=zms(ip3)+zn
 822  continue
        xyzmax=xpetit/xzprec
        do ip = 1,nbpts
         if (ims(ip).ne.0) then
          xyzmax=max(xyzmax,abs(xms(ip)))
          xyzmax=max(xyzmax,abs(yms(ip)))
          xyzmax=max(xyzmax,abs(zms(ip)))
         endif
        enddo
        do 821 ip = 1,nbpts
         if(ims(ip).ne.0) then
          sn = xms(ip)*xms(ip)+yms(ip)*yms(ip)+zms(ip)*zms(ip)
          sn = sqrt(abs(sn))
          if (sn.lt.xyzmax*1d-10) then
            interr(1)=ip
            call erreur(1151)
            return
          endif
         xms(ip)=xms(ip)/sn
         yms(ip)=yms(ip)/sn
         zms(ip)=zms(ip)/sn
**       write(6,*) 'ip ',ip,xms(ip),yms(ip),zms(ip)
         endif
 821  continue
        if(ipt1.itypel.ne.22) call erreur(16)
**      write(6,*) 'ipt8.itypel ',ipt8.itypel, ipt7
        if(ierr.ne.0) return
*
*  Nombre d'iterations pour la detection de la direction du contact
      nbiter=16
*
* Relation du type 3 : noeud-triangle
*-------------------------------------
*
*  creation du meleme associe a la relation
*
      nbsous=0
      nbref =0
      nbnn  =5
      nbelem=incnel
      segini meleme
      itypel=22
      irigel(1,nrigel)=meleme
*
*  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 10 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)
 10   continue
      segdes,descr
      irigel(3,nrigel) = descr
*
*  creation du segment xmatri
*
      nelrig = incnel
      segini,xmatri
      irigel(4,nrigel) = xmatri
*
*  ce qu'on cree est unilateral
*
      irigel(6,nrigel)=1
*
*  ce qu'on cree est symetrique
*
      irigel(7,nrigel)=0
*
*  Nombre d'elements dans la rigidite du type 3
      nelri3 = 0
*
*  Preparation du zonage de ipt1
*
      nbzx=(xmax-xmin)/tamax
      nbzy=(ymax-ymin)/tamax
      nbzz=(zmax-zmin)/tamax
*
*  preconditionnement de ipt1: mfopa4
*
      segini mfopa4
      xmult=3.1415926*(xmax-xmin)
      ymult=2.7182818*(ymax-ymin)
      zmult=1.*(zmax-zmin)
      tmult=sqrt(xmult**2+ymult**2+zmult**2)
      xmult=xmult/tmult
      ymult=ymult/tmult
      zmult=zmult/tmult
      tmin=xgrand
      tmax=-xgrand
      do  j=1,nbel1
       ip=ipt1.num(2,j)
        ipv = (ip-1)*idimp1
        xp = xcoor(ipv+1)
        yp = xcoor(ipv+2)
        zp = xcoor(ipv+3)
        tco(j)=xp*xmult+yp*ymult+zp*zmult
        tmin=min(tco(j),tmin)
        tmax=max(tco(j),tmax)
        ipnum(j)=ip
        ico(j)=j
      enddo
      nbzt=(tmax-tmin)/tamax
      nbzt=max(nbzt,1)
      nbzt=min(nbel1,nbzt)
**    write(6,*) ' nbzt ',nbzt
 
*
*  trier selon x
*
      call triflo(tco,xw,ico,iw,nbel1)
*
*  indexer
*     
      indt=nbzt+1
      segini mfopa5
      do i=nbel1,1,-1
        id=nbzt*(tco(i)-tmin)/(tmax-tmin)+1
        ind(id)=i
      enddo
      do i=1,nbzt
        if (ind(i+1).eq.0) ind(i+1)=ind(i)
        if(ind(i+1).lt.ind(i)) call erreur(5)
      enddo
 
 
*  Boucle sur le maillage support du contact(frottement)
*
**    write(6,*) ' boucle 19 nbel6 ',nbel6
      DO 19 iel6=1,nbel6
*
        xjr = 0d0
        ip1=ipt6.num(1,iel6)
        ip2=ipt6.num(2,iel6)
        ip3=ipt6.num(3,iel6)
*  Recuperation des coordonees des noeuds du triangle
        ipv = (ip1-1)*idimp1
        xp1 = xcoor(ipv+1)
        yp1 = xcoor(ipv+2)
        zp1 = xcoor(ipv+3)
        ipv = (ip2-1)*idimp1
        xp2 = xcoor(ipv+1)
        yp2 = xcoor(ipv+2)
        zp2 = xcoor(ipv+3)
        ipv = (ip3-1)*idimp1
        xp3 = xcoor(ipv+1)
        yp3 = xcoor(ipv+2)
        zp3 = xcoor(ipv+3)
*  Centre de gravite du triangle
        xg = (xp1+xp2+xp3) /3.d0
        yg = (yp1+yp2+yp3) /3.d0
        zg = (zp1+zp2+zp3) /3.d0
* critere de distance
        d1 = (xg-xp1)**2 + (yg-yp1)**2 + (zg-zp1)**2
        d2 = (xg-xp2)**2 + (yg-yp2)**2 + (zg-zp2)**2
        d3 = (xg-xp3)**2 + (yg-yp3)**2 + (zg-zp3)**2
*  Triangle un peu plus grand pour les tests
**      scale=1.00 + xszpre
        scale=1.00 + 1D-4
        dist2  = max(d1,d2,d3)
        dist = sqrt(dist2)*scale
        sqdist4 = dist * 4  
        xp1e=xg+(xp1-xg)*scale
        yp1e=yg+(yp1-yg)*scale
        zp1e=zg+(zp1-zg)*scale
        xp2e=xg+(xp2-xg)*scale
        yp2e=yg+(yp2-yg)*scale
        zp2e=zg+(zp2-zg)*scale
        xp3e=xg+(xp3-xg)*scale
        yp3e=yg+(yp3-yg)*scale
        zp3e=zg+(zp3-zg)*scale
  25    continue
*   rechercher le noeud a tester dans le deuxieme maillage qui est sous forme mult avec en
*   deuxieme position le point physique
**      write(6,*) ' boucle 20 ipt1 ',ipt1.num(/2)
        xgm = xg -sqdist4
        xgp = xg +sqdist4
        ygm = yg -sqdist4
        ygp = yg +sqdist4
        zgm = zg -sqdist4
        zgp = zg +sqdist4
*
*  calcul zone du centre de gravite
*
        tc=xg*xmult+yg*ymult+zg*zmult
*
        xmn = (xms(ip1)+xms(ip2)+xms(ip3))*X1s3
        ymn = (yms(ip1)+yms(ip2)+yms(ip3))*X1s3
        zmn = (zms(ip1)+zms(ip2)+zms(ip3))*X1s3
        xzonag = sqdist4*sqr3
        tjeup=0.d0
        tjeum=0.d0
 
        if (MELVA1.ne.0) then
*  xmult ymult et zmult doivent etre positifs
          tjeum=xmult*(xjmin*xmn)+ymult*(xjmin*ymn)+zmult*(xjmin*zmn)
          tjeup=xmult*(xjmax*xmn)+ymult*(xjmax*ymn)+zmult*(xjmax*zmn)
        endif
*       
        izg=nbzt*((tc-xzonag-tjeum)-tmin)/(tmax-tmin)+1
        izg=max(izg,1)
        izg=min(izg,indt)
* debut de zone
        indb=ind(izg)
        do 20 iz=indb,nbel1
          if(tco(iz).lt.(tc-tjeum-xzonag)) goto 20
          if(tco(iz).gt.(tc-tjeup+xzonag)) then
**          write(6,*) ' sortie pour ',iz, ' en ',iz-indb+1
            goto 18
          endif
        iel1 = ico(iz)
        jp = ipnum(iel1)
*  Verification que pas relation du point sur lui meme
        if (jp.eq.ip1) goto 20
        if (jp.eq.ip2) goto 20
        if (jp.eq.ip3) goto 20
*  verification rapide en norme L1 d'abord
        ipv = (jp-1)*idimp1
        xp = xcoor(ipv+1)
        yp = xcoor(ipv+2)
        zp = xcoor(ipv+3)
*
        xpp=xp
        ypp=yp
        zpp=zp
        if (MELVA1.ne.0) then
          nel1 = min (iel1,nelj)
          xjr = melva1.velche(nptelj,nel1)
          xpp=xpp-xjr*xmn
          ypp=ypp-xjr*ymn
          zpp=zpp-xjr*zmn
        endif
        if(xpp.lt.xgm.or.xpp.gt.xgp) goto 20
        if(ypp.lt.ygm.or.ypp.gt.ygp) goto 20
        if(zpp.lt.zgm.or.zpp.gt.zgp) goto 20
*
*  Verification si autre point dans la zone de selection
*  verif distance au centre de gravite
        dp = ((xg-xpp)**2 + (yg-ypp)**2 + (zg-zpp)**2)
**      write(6,*) 'dp dist2',dp,dist2,xjr
***        if (dp.gt.xzonag**2) then
***         goto 20
***        endif
C*DBG        write(ioimp,*) 'contact test distance ok',dp,d1,d2,d3
 
*  verif position par rapport aux cotes
 
*  cote 1 2
        x12 = xp2 - xp1
        y12 = yp2 - yp1
        z12 = zp2 - zp1
        sv = sqrt(x12**2+y12**2+z12**2)
        xv=x12/sv
        yv=y12/sv
        zv=z12/sv
*  normale locale (1 et 2)
        xnl=xms(ip1)+xms(ip2)
        ynl=yms(ip1)+yms(ip2)
        znl=zms(ip1)+zms(ip2)
 
*  vecteur reference
        xn=y12*znl-z12*ynl
        yn=z12*xnl-x12*znl
        zn=x12*ynl-y12*xnl
        dn  = sqrt(xn*xn+yn*yn+zn*zn)
        scal = (xpp-xp1e)*xv + (ypp-yp1e)*yv + (zpp-zp1e)*zv
        xm = xp1e + scal*xv
        ym = yp1e + scal*yv
        zm = zp1e + scal*zv
        dpm = sqrt(abs((xpp-xm)**2+(ypp-ym)**2+(zpp-zm)**2))
        scal = (xpp-xm)*xn + (ypp-ym)*yn + (zpp-zm)*zn
**      if (dpm.lt.dist*xszpre) write(6,*) ' pt sur 1 2',dpm
        if (dpm.gt.dist*xszpre.and.scal.gt.0.707d0*dn*dpm) then
*        write(6,*) ' 1 dpm scal ',dpm,scal
        goto 20
        endif
        dpm=max(xpetit,dpm)
*       write(ioimp,*) 'contact test position 1 ok',
*    &                 scal/(dn*dpm),scal,dn,dpm
*
*  cote 2 3
        x23 = xp3 - xp2
        y23 = yp3 - yp2
        z23 = zp3 - zp2
        sv = sqrt(x23**2+y23**2+z23**2)
        xv=x23/sv
        yv=y23/sv
        zv=z23/sv
*  normale locale (2 et 3)
        xnl=xms(ip2)+xms(ip3)
        ynl=yms(ip2)+yms(ip3)
        znl=zms(ip2)+zms(ip3)
 
*  vecteur reference
        xn=y23*znl-z23*ynl
        yn=z23*xnl-x23*znl
        zn=x23*ynl-y23*xnl
        dn  = sqrt(xn*xn+yn*yn+zn*zn)
        scal = (xpp-xp2e)*xv + (ypp-yp2e)*yv + (zpp-zp2e)*zv
        xm = xp2e + scal*xv
        ym = yp2e + scal*yv
        zm = zp2e + scal*zv
        dpm = sqrt(abs((xpp-xm)**2+(ypp-ym)**2+(zpp-zm)**2))
        scal = (xpp-xm)*xn + (ypp-ym)*yn + (zpp-zm)*zn
**      if (dpm.lt.dist*xszpre) write(6,*) ' pt sur 2 3',dpm
        if (dpm.gt.dist*xszpre.and.scal.gt.0.707d0*dn*dpm) then
*        write(6,*) ' 2 dpm scal ',dpm,scal
        goto 20
        endif
        dpm=max(xpetit,dpm)
*       write(ioimp,*) 'contact test position 2 ok',
*    &                 scal/(dn*dpm),scal,dn,dpm
*
*  cote 3 1
        x31 = xp1 - xp3
        y31 = yp1 - yp3
        z31 = zp1 - zp3
        sv = sqrt(x31**2+y31**2+z31**2)
        xv=x31/sv
        yv=y31/sv
        zv=z31/sv
*  normale locale (3 et 1)
        xnl=xms(ip3)+xms(ip1)
        ynl=yms(ip3)+yms(ip1)
        znl=zms(ip3)+zms(ip1)
*  vecteur reference
        xn=y31*znl-z31*ynl
        yn=z31*xnl-x31*znl
        zn=x31*ynl-y31*xnl
        dn  = sqrt(xn*xn+yn*yn+zn*zn)
        scal = (xpp-xp3e)*xv + (ypp-yp3e)*yv + (zpp-zp3e)*zv
        xm = xp3e + scal*xv
        ym = yp3e + scal*yv
        zm = zp3e + scal*zv
        dpm = sqrt(abs((xpp-xm)**2+(ypp-ym)**2+(zpp-zm)**2))
        scal = (xpp-xm)*xn + (ypp-ym)*yn + (zpp-zm)*zn
**      if (dpm.lt.dist*xszpre) write(6,*) ' pt sur 3 1',dpm
        if (dpm.gt.dist*xszpre.and.scal.gt.0.707d0*dn*dpm) then
*        write(6,*) ' 3 dpm scal ',dpm,scal
        goto 20
        endif
        dpm=max(xpetit,dpm)
*       write(ioimp,*) 'contact test position 2 ok',
*    &                 scal/(dn*dpm),scal,dn,dpm
*
* on a un point ou imposer la relation
C*DBG       write(ioimp,*) 'contact potentiel '
* Quelle est la relation ???
*
*   direction de la relation = normale au plan (123)
*
*  normale reelle
        xnr = y12*z23 - z12*y23
        ynr = z12*x23 - x12*z23
        znr = x12*y23 - y12*x23
        dnr = sqrt(xnr*xnr+ynr*ynr+znr*znr)
        if (dnr.lt.dist*xszpre.AND.iimpi.ne.0) write(ioimp,*) ' pb dnr'
        xnr = xnr / dnr
        ynr = ynr / dnr
        znr = znr / dnr
* normale ponderee
        xn=xms(ip1)+xms(ip2)+xms(ip3)
        yn=yms(ip1)+yms(ip2)+yms(ip3)
        zn=zms(ip1)+zms(ip2)+zms(ip3)
        dn = sqrt(xn*xn+yn*yn+zn*zn)
        if (.not.(dn.lt.1d-10).and. .not.(dn.gt.1d-10).AND.iimpi.ne.0)
     &    write(ioimp,*) 'Prob 4.1 - impo32'
        if (abs(dn).le.xpetit.AND.iimpi.ne.0)
     &    write(ioimp,*) 'Prob 4.2 - impo32'
        xn = xn / dn
        yn = yn / dn
        zn = zn / dn
*
*   Distance (jeu) du point jp au plan de contact
*   Puis calcul de la projection sur le plan de contact
*
**      write(6,*) 'ip1 xms ',ip1,xms(ip1),yms(ip1),zms(ip1)
**      write(6,*) 'ip2 xms ',ip2,xms(ip2),yms(ip2),zms(ip2)
**      write(6,*) 'ip3 xms ',ip3,xms(ip3),yms(ip3),zms(ip3)
**      write(6,*) 'xn      ',xn,yn,zn
        iter = 0
        xjeu  = (xp-xg)*xnr + (yp-yg)*ynr + (zp-zg)*znr
  21    continue
        angn = xn * xnr + yn*ynr + zn*znr
        if (angn.le.xpetit.AND.iimpi.ne.0)
     &    write(ioimp,*) 'angn negatif ',angn
        if(angn.le.xpetit) goto 20
        xjeuv = xjeu  / angn
*
*   Recherche de ses coordonnées barycentriques
*   qui sont les surfaces des sous triangles
*
        xq = xp - xjeuv*xn
        yq = yp - xjeuv*yn
        zq = zp - xjeuv*zn
        xb1 = (yq-yp2)*(zq-zp3)-(zq-zp2)*(yq-yp3)
        yb1 = (zq-zp2)*(xq-xp3)-(xq-xp2)*(zq-zp3)
        zb1 = (xq-xp2)*(yq-yp3)-(yq-yp2)*(xq-xp3)
        xb2 = (yq-yp3)*(zq-zp1)-(zq-zp3)*(yq-yp1)
        yb2 = (zq-zp3)*(xq-xp1)-(xq-xp3)*(zq-zp1)
        zb2 = (xq-xp3)*(yq-yp1)-(yq-yp3)*(xq-xp1)
        xb3 = (yq-yp1)*(zq-zp2)-(zq-zp1)*(yq-yp2)
        yb3 = (zq-zp1)*(xq-xp2)-(xq-xp1)*(zq-zp2)
        zb3 = (xq-xp1)*(yq-yp2)-(yq-yp1)*(xq-xp2)
        b1  = xb1*xnr + yb1*ynr + zb1*znr
        b2  = xb2*xnr + yb2*ynr + zb2*znr
        b3  = xb3*xnr + yb3*ynr + zb3*znr
        bt  = b1 + b2 + b3
*  normalement pas utile
*  element retourne a cause des grands deplacements
        if (bt.lt.xpetit) then
        write (ioimp,*) ' bt negatif dans impo32 '
           call soucis(719)
**         bt=xpetit*2.d0
           goto 20
        endif
        if (bt.lt.0.d0.AND.iimpi.ne.0) then
        write (ioimp,*) ' bt negatif dans impo32 '
         print *,'xp1 yp1 zp1',xp1,yp1,zp1
         print *,'xp2 yp2 zp2',xp2,yp2,zp2
         print *,'xp3 yp3 zp3',xp3,yp3,zp3
         print *,'xp yp zp',xp,yp,zp
         print *,'xn yn zn',xn,yn,zn
         print *,'b1 b2 b3 bt',b1,b2,b3,bt
*        goto 20
        endif
        bsurf = bt / 2
        if (.not.(bt.lt.1d-10) .and. .not.(bt.gt.1d-10).AND.iimpi.ne.0)
     &    write(ioimp,*) 'Prob 5.1 - impo32'
        if (abs(bt).le.xpetit) write(ioimp,*) 'Prob 5.2 - impo32'
        b1 = b1 / bt
        b2 = b2 / bt
        b3 = b3 / bt
*  recalcul de la direction et du jeu en fonction des normales aux sommets
*  on arete l'interpolation de la normale au bord de l'element
        bb1=max(b1,0.d0)
        bb2=max(b2,0.d0)
        bb3=max(b3,0.d0)
        bm = bb1 + bb2 + bb3
        bb1 = bb1 / bm
        bb2 = bb2 / bm
        bb3 = bb3 / bm
        xnp = xn
        ynp = yn
        znp = zn
        xn =   bb1*xms(ip1)+ bb2*xms(ip2)+ bb3*xms(ip3)
        yn =   bb1*yms(ip1)+ bb2*yms(ip2)+ bb3*yms(ip3)
        zn =   bb1*zms(ip1)+ bb2*zms(ip2)+ bb3*zms(ip3)
        sn = sqrt (xn*xn + yn*yn + zn*zn)
        xn=xn/sn
        yn=yn/sn
        zn=zn/sn
        diff=abs(xn-xnp)+abs(yn-ynp)+abs(zn-znp)
*       write (6,*) ' iter b* *n ',iter,b1,b2,b3,xn,yn,zn
* recalcul en fonction de la nouvelle normale
        iter=iter+1
        if (iter.gt.64) then
**        write(6,*) ' impo32 diff ',diff
        goto 20
        endif
        if (diff.gt.1d-10) goto 21
**      write(6,*) ' impo32 iter ',iter
**        write (6,*) ' iter b* *n ',iter,b1,b2,b3,xn,yn,zn
 
 
* si on a deja traverse, les trois coordonnees barycentriques doivent etre positives
**      if (xjeuv-xjr.lt.-xszpre*dist) then
**        if (b1.lt.-xszpre) goto 20
**        if (b2.lt.-xszpre) goto 20
**        if (b3.lt.-xszpre) goto 20
**      endif
*  Si on est en dehors, on projette sur l'arete (ou pas)
**      if (b1.lt.0.d0.or.b2.lt.0.d0.or.b3.lt.0.d0) then
**        xq=xp1*bb1+xp2*bb2+xp3*bb3
**        yq=yp1*bb1+yp2*bb2+yp3*bb3
**        zq=zp1*bb1+zp2*bb2+zp3*bb3
**        xnn=xp-xq
**        ynn=yp-yq
**        znn=zp-zq
**        snn=sqrt(abs(xnn*xnn+ynn*ynn+znn*znn))
**        if (snn.lt.xszpre*dist) write (6,*) ' snn petit ',snn
*  sinon on prend la direction reelle
**        if (snn.gt.xszpre*dist) then
**           xn=xnn/sn
**           yn=ynn/sn
**           zn=znn/sn
**        endif
**      endif
        xjeu1 = (xp-xp1)*xn + (yp-yp1)*yn + (zp-zp1)*zn
        xjeu2 = (xp-xp2)*xn + (yp-yp2)*yn + (zp-zp2)*zn
        xjeu3 = (xp-xp3)*xn + (yp-yp3)*yn + (zp-zp3)*zn
*****pv        xjeuv = xjeu1 * bb1 + xjeu2 * bb2 + xjeu3 * bb3
**      write (6,*) 'b1 b2 b3',b1,b2,b3
**      write (6,*) 'xjeu xjeuv ',xjeu,xjeuv,xjeu1,xjeu2,xjeu3
        xjeu = xjeuv - xjr
**      write (6,*) ' normale ',xn,yn,zn,' jeu ',xjeu1,xjeu2,xjeu3
*  verif bon cote (a la tolerance pres)
**      write (6,*) ' xjeu ',xjeu,dist,iel6
**        if (xjeu.lt.-1*dist) goto 20
**        if (xjeu.gt.3*dist) goto 20
*  verif compatibilite avec la normale au poin impactant
      if (xms(jp)*xn+yms(jp)*yn+zms(jp)*zn.gt. -0.0d0) then
*         write (6,*) ' impo32 normales incompatibes 1 ',
*    >  jp,xjeu,dpm,dist
       goto 20
      endif
 
 
 
C*DBG        write(ioimp,*) ' b1 b2 b3 ',b1,b2,b3
 1954   continue
* points a l'exterieur?
* on met une ponderation et un rayon d'acceptation
        pond =  (1D4+1) - 1D4 * bm
        if (pond.le.0.d0) goto 20
        pond=max(pond,0.d0)
        pond=min(pond,1.d0)
 
**       if (xjeu.lt.0.d0)
**   >    write (6,*) 'pt traverse',jp,b1,b2,b3,xjeuv-xjr,dist
 
*
*  Ajout d'une relation noeud-triangle
*
        nelri3 = nelri3 + 1
        nelch = nelch + 1
*
*   on ajuste les differents segments si necesaire
*
        if (nelri3.gt.nelrig) then
          nelrig = nelrig + incnel
          segadj,xmatri
          nbelem = nbelem + incnel
          nbnn = 5
          segadj,meleme
          nbnn = 1
          segadj,ipt7
          n = n + incnel
          segadj,mpoval
        endif
*
*  Mise a jour du meleme
**      write(6,*) 'ipt8.num(/2)',ipt8.num(/2)
        num(1,nelri3) = ipt1.num(1,iel1)
        num(2,nelri3) = ip1
        num(3,nelri3) = ip2
        num(4,nelri3) = ip3
        num(5,nelri3) = jp
**      write(6,*) ' nouveau meleme',num(1,nelri3),num(2,nelri3),
**   >  num(3,nelri3),num(4,nelri3),num(5,nelri3)                                                  *
*  Mise a jour de xmatri
*  lambda
        re( 1,1,nelri3) = 0d0
*  ip1
        re( 2,1,nelri3) = xn * bb1 * bsurf * pond
        re( 3,1,nelri3) = yn * bb1 * bsurf * pond
        re( 4,1,nelri3) = zn * bb1 * bsurf * pond
*  ip2
        re( 5,1,nelri3) = xn * bb2 * bsurf * pond
        re( 6,1,nelri3) = yn * bb2 * bsurf * pond
        re( 7,1,nelri3) = zn * bb2 * bsurf * pond
*  ip3
        re( 8,1,nelri3) = xn * bb3 * bsurf * pond
        re( 9,1,nelri3) = yn * bb3 * bsurf * pond
        re(10,1,nelri3) = zn * bb3 * bsurf * pond
*  jp
        re(11,1,nelri3) = -xn * bsurf * pond
        re(12,1,nelri3) = -yn * bsurf * pond
        re(13,1,nelri3) = -zn * bsurf * pond
*  on transpose
        do 30 ic = 2, nligrp
          re(1,ic,nelri3) = re(ic,1,nelri3)
 30     continue
*  le reste est nul
*
*  remplissage du champoint de depi (jeu)
*
        ipt7.num(1,nelch) = ipt1.num(1,iel1)
**      write(6,*) ' ipt1.num(1,iel1) ',ipt1.num(1,iel1)
        vpocha(nelch,1)   = xjeu  * bsurf * pond
        vpocha(nelch,2)   = bsurf * pond
        IF (MELVA2.NE.0) THEN
          NEL2 = min (iel1,NELA)
          VPOCHA(nelch,3) = max(MELVA2.VELCHE(NPTELA,NEL2),0.D0)*
     >      bsurf * pond
        ENDIF
**        write(ioimp,*) ' jeu ',xjeu,ip1,ip2,ip3,jp,b1,b2,b3,xn,yn,zn
*
 20   CONTINUE
 18   CONTINUE
 19   CONTINUE
*
* Ajustement au plus juste puis desactivation des segments lies
* a la rigidite du type 3
      if (nelri3.ne.nelrig) then
        nelrig = nelri3
        segadj,xmatri
        nbelem = nelri3
        nbnn = 5
        segadj,meleme
      endif
      segdes,xmatri
*
**         write(ioimp,*) ' nb relation type 3 ',nelri3
**         write(6,fmt='(10i8)') (num(5,nelr),nelr=1,nelri3)
 
 
 
 
 
C*DBG      if (nelri3.eq.0) irigel(6,nrigel)=0
*
*  Fin du traitement de la formulation standard/forte
* ----------------------------------------------------
 300  CONTINUE
* Destruction des segments locaux devenus inutiles
      segsup mfopa3,mfopa4,mfopa5
 
      GOTO 1000
 
*=======================================================================
*= Formulation "faible" 3D : relation triangle-triangle
*=======================================================================
* Element de contact 3D a 7 noeuds (1+6)
 500  CONTINUE
*
*  Petit segment de travail
      segini,mfaible
*
* Relation du type 0 : triangle-triangle (faible)
*---------------------------------------
*
*  Creation du meleme associe a la relation
*
      nbnn   = 7
      nbelem = incnel
      nbsous = 0
      nbref  = 0
      segini,meleme
      itypel = 22
      irigel(1,nrigel) = meleme
*
*  Creation du descriptif commun a toutes les raideurs
*
      nligrp = 19
      nligrd = nligrp
      segini,descr
      lisinc(1) = 'LX  '
      lisdua(1) = 'FLX '
      noelep(1) = 1
      noeled(1) = 1
      do 510 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+5)/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)
 510  continue
      segdes,descr
      irigel(3,nrigel) = descr
*
*     creation du segment xmatri
*
      nelrig = incnel
      segini,xmatri
      irigel(4,nrigel) = xmatri
*
*     ce qu'on cree est unilateral
*
      irigel(6,nrigel) = 1
*
*     ce qu'on cree est symetrique
*
      irigel(7,nrigel) = 0
*
*  Nombre d'elements crees dans meleme=irigel(1,nrigel), ipt7 et mpoval
      nelri0 = 0
*
*  Boucle sur les elements du maillage de contact/frottement "faible"
*
      DO 519 iel8=1,nbel8
*
        xjr = 0d0
        ip1 = ipt8.num(1,iel8)
        ip2 = ipt8.num(2,iel8)
        ip3 = ipt8.num(3,iel8)
* Definition du triangle T1(ip1,ip2,ip3)
* - Coordonnees des noeuds de T1
        ipv = (ip1-1)*idimp1
        xp1 = xcoor(ipv+1)
        yp1 = xcoor(ipv+2)
        zp1 = xcoor(ipv+3)
        ipv = (ip2-1)*idimp1
        xp2 = xcoor(ipv+1)
        yp2 = xcoor(ipv+2)
        zp2 = xcoor(ipv+3)
        ipv = (ip3-1)*idimp1
        xp3 = xcoor(ipv+1)
        yp3 = xcoor(ipv+2)
        zp3 = xcoor(ipv+3)
* - Barycentre de T1
        xgT1 = (xp1 + xp2 + xp3) * X1s3
        ygT1 = (yp1 + yp2 + yp3) * X1s3
        zgT1 = (zp1 + zp2 + zp3) * X1s3
* - Normale au triangle T1
        xnT1 = (yp1-yp2)*(zp2-zp3) - (zp1-zp2)*(yp2-yp3)
        ynT1 = (zp1-zp2)*(xp2-xp3) - (xp1-xp2)*(zp2-zp3)
        znT1 = (xp1-xp2)*(yp2-yp3) - (yp1-yp2)*(xp2-xp3)
        dnT1 = SQRT(xnT1*xnT1+ynT1*ynT1+znT1*znT1)
        IF (.not.(dnT1.lt.1d-10) .and. .not.(dnT1.gt.1d-10).AND.
     &      iimpi.ne.0) write(ioimp,*) 'FAIBle - 1.1 - impo32'
        IF (abs(dnT1).le.xpetit) write(ioimp,*) 'FAIBle - 1.2 - impo32'
        xnT1 = xnT1 / dnT1
        ynT1 = ynT1 / dnT1
        znT1 = znT1 / dnT1
C*DBG        write(ioimp,*) 'Normale au plan T1',xnT1,ynT1,znT1,dnT1
* - Rayon de la sphere "englobante" de l'element (centre=barycentre)
        d1 = (xgT1-xp1)**2 + (ygT1-yp1)**2 + (zgT1-zp1)**2
        d2 = (xgT1-xp2)**2 + (ygT1-yp2)**2 + (zgT1-zp2)**2
        d3 = (xgT1-xp3)**2 + (ygT1-yp3)**2 + (zgT1-zp3)**2
        RayT1 = SQRT(max(d1,d2,d3))
      DO 520 iel6 = 1, nbel6
        if (MELVA1.ne.0) then
           nel1 = min (iel6,nelj)
           xjr = melva1.velche(nptelj,nel1)
        endif
        jp1 = ipt6.num(1,iel6)
* - Les noeuds 2 et 3 sont intervertis (fait aupravant dans impo31)
*   Cela permet aux triangles T1 et T2 d etre orientes dans le meme sens.
        jp3 = ipt6.num(2,iel6)
        jp2 = ipt6.num(3,iel6)
C*DBG        write(ioimp,*) iel,ip1,ip2,ip3,jp1,jp2,jp3,ipt1.num(1,iel)
*
* Verification (provisoire) qu'il n'y a pas de noeud double
        if (ip1.eq.jp1) goto 520
        if (ip1.eq.jp2) goto 520
        if (ip1.eq.jp3) goto 520
        if (ip2.eq.jp1) goto 520
        if (ip2.eq.jp2) goto 520
        if (ip2.eq.jp3) goto 520
        if (ip3.eq.jp1) goto 520
        if (ip3.eq.jp2) goto 520
        if (ip3.eq.jp3) goto 520
*
*
* Definition du triangle T2(jp1,jp2,jp3)
* - Coordonnees des noeuds de T2
        ipv = (jp1-1)*idimp1
        xq1 = xcoor(ipv+1)
        yq1 = xcoor(ipv+2)
        zq1 = xcoor(ipv+3)
        ipv = (jp2-1)*idimp1
        xq2 = xcoor(ipv+1)
        yq2 = xcoor(ipv+2)
        zq2 = xcoor(ipv+3)
        ipv = (jp3-1)*idimp1
        xq3 = xcoor(ipv+1)
        yq3 = xcoor(ipv+2)
        zq3 = xcoor(ipv+3)
* - Barycentre de T2
        xgT2 = (xq1 + xq2 + xq3) * X1s3
        ygT2 = (yq1 + yq2 + yq3) * X1s3
        zgT2 = (zq1 + zq2 + zq3) * X1s3
* - Normale au triangle T2
        xnT2 = (yq1-yq2)*(zq2-zq3) - (zq1-zq2)*(yq2-yq3)
        ynT2 = (zq1-zq2)*(xq2-xq3) - (xq1-xq2)*(zq2-zq3)
        znT2 = (xq1-xq2)*(yq2-yq3) - (yq1-yq2)*(xq2-xq3)
        dnT2 = SQRT(xnT2*xnT2+ynT2*ynT2+znT2*znT2)
        IF (.not.(dnT2.lt.1d-10) .and. .not.(dnT2.gt.1d-10).AND.
     &      iimpi.ne.0) write(ioimp,*) 'FAIBle - 2.1 - impo32'
        IF (abs(dnT2).le.xpetit.AND.iimpi.ne.0)
     &    write(ioimp,*) 'FAIBle - 2.2 - impo32'
        xnT2 = xnT2 / dnT2
        ynT2 = ynT2 / dnT2
        znT2 = znT2 / dnT2
C*DBG        write(ioimp,*) 'Normale au plan T2',xnT2,ynT2,znT2,dnT2
* - Rayon de la sphere "englobante" de l'element (centre=barycentre)
        d1 = (xgT2-xq1)**2 + (ygT2-yq1)**2 + (zgT2-zq1)**2
        d2 = (xgT2-xq2)**2 + (ygT2-yq2)**2 + (zgT2-zq2)**2
        d3 = (xgT2-xq3)**2 + (ygT2-yq3)**2 + (zgT2-zq3)**2
        RayT2 = SQRT(max(d1,d2,d3))
 
* Orientation respective correcte des 2 normales
*  Les triangles T1 et T2 etant decrits dans le meme sens ("prisme"),
*  le produit scalaire de leur normale doit etre positif ou nul !
        scal = xnT1 * xnT2 + ynT1 * ynT2 + znT1 * znT2
C*DBG        write(ioimp,*) 'Prod.scal des normales',scal
        if (scal.lt.0.) goto 520
*
* Distance entre les centres de gravite
* Les spheres "englobantes" doivent s'intersecter pour contact potentiel
        dist = SQRT((xgT2-xgT1)**2 + (ygT2-ygT1)**2 + (zgT2-zgT1)**2)
*       IF (dist*4.GT.RayT1+RayT2) goto 520
*
* Definition du plan de contact (surface intermediaire)
* - Normale (unitaire) du plan de contact
        xnC = xnT1 + xnT2
        ynC = ynT1 + ynT2
        znC = znT1 + znT2
        dnC = SQRT(xnC*xnC + ynC*ynC + znC*znC)
        IF (.not.(dnC.lt.1d-10).and. .not.(dnC.GT.1D-10).AND.iimpi.ne.0)
     &    write(ioimp,*) 'FAIBle - 3.1 - impo32'
        IF (abs(dnC).le.xpetit.AND.iimpi.ne.0)
     &    write(ioimp,*) 'FAIBle - 3.2 - impo32'
        xnC = xnC / dnC
        ynC = ynC / dnC
        znC = znC / dnC
C*DBG        write(ioimp,*) 'Normale au plan Contact',xnC,ynC,znC,dnC
* 
* Criteres d'acceptation de l'element de contact
        RayTT = RayT1 + RayT2
* - Les barycentres des triangles T1 et T2 sont-ils suffisamment proches :
*       distance suivant la direction de contact (prise en compte du jeu)
        VXB12 = xgT2 - xgT1 
        VYB12 = ygT2 - ygT1 
        VZB12 = zgT2 - zgT1 
        distn = (VXB12*xnC)+(VYB12*ynC)+(VZB12*znC)
        test1 = distn
        if (MELVA1.NE.0) then
         test1 = test1 - xjr
        endif
        if (test1.GT.RayTT) GOTO 520
*
*       distance dans un plan dont la normale est la direction de contact
        distt = SQRT(dist**2 - distn**2)
        if (distt.GT.RayTT) GOTO 520
 
* - Point definissant le plan de contact
*   = milieu des barycentres des triangles T1 et T2
        xgC = (xgT1 + xgT2) * 0.5
        ygC = (ygT1 + ygT2) * 0.5
        zgC = (zgT1 + zgT2) * 0.5
* - "Distance" de reference au plan de contact
        dgC = xgC*xnC + ygC*ynC + zgC*znC
*
* Distance "signee" des triangles T1 et T2 au plan de contact
        dp1C = (xp1*xnC + yp1*ynC + zp1*znC) - dgC
        dp2C = (xp2*xnC + yp2*ynC + zp2*znC) - dgC
        dp3C = (xp3*xnC + yp3*ynC + zp3*znC) - dgC
        dq1C = (xq1*xnC + yq1*ynC + zq1*znC) - dgC
        dq2C = (xq2*xnC + yq2*ynC + zq2*znC) - dgC
        dq3C = (xq3*xnC + yq3*ynC + zq3*znC) - dgC
* Projection des triangles T1 et T2 sur le plan de contact
* T1C = triangle T1 projete sur plan de Contact
        cPT1C(1,1) = xp1 - dp1C * xnC
        cPT1C(1,2) = yp1 - dp1C * ynC
        cPT1C(1,3) = zp1 - dp1C * znC
        cPT1C(2,1) = xp2 - dp2C * xnC
        cPT1C(2,2) = yp2 - dp2C * ynC
        cPT1C(2,3) = zp2 - dp2C * znC
        cPT1C(3,1) = xp3 - dp3C * xnC
        cPT1C(3,2) = yp3 - dp3C * ynC
        cPT1C(3,3) = zp3 - dp3C * znC
* T2C = triangle T2 projete sur plan de Contact
        cPT2C(1,1) = xq1 - dq1C * xnC
        cPT2C(1,2) = yq1 - dq1C * ynC
        cPT2C(1,3) = zq1 - dq1C * znC
        cPT2C(2,1) = xq2 - dq2C * xnC
        cPT2C(2,2) = yq2 - dq2C * ynC
        cPT2C(2,3) = zq2 - dq2C * znC
        cPT2C(3,1) = xq3 - dq3C * xnC
        cPT2C(3,2) = yq3 - dq3C * ynC
        cPT2C(3,3) = zq3 - dq3C * znC
*
* On determine quelle est la composante maximale de la normale pour
* projeter les triangles T1C et T2C sur un plan parallele aux axes
* ("le plus orthogonal" a cette normale), ce qui maximise leur aire.
        r_x = abs(xnC)
        r_y = abs(ynC)
        r_z = abs(znC)
        if (r_x.gt.r_y) then
          inC = 1
          if (r_x.lt.r_z) inC = 3
        else
          inC = 2
          if (r_y.lt.r_z) inC = 3
        endif
* Projection des triangles T1C et T2C sur ce plan // aux axes (A)
* T1A = triangle T1C projete sur ce plan (A)
* T2A = triangle T2C projete sur ce plan (A)
* On prend soin d'orienter les triangles T1A et T2A dans le sens
* direct (en tenant compte du signe de la composante normale !)
        GOTO (531,532,533), inC
        call erreur(5)
        return
* Normale NC selon "x" : A =(y,z)
 531    if (xnC.lt.0.) inC = -inC
        inX = 2
        inY = 3
        goto 534
* Normale NC selon "y" : A =(x,z)
 532    if (ynC.lt.0.) inC = -inC
        inX = 1
        inY = 3
        goto 534
* Normale NC selon "z" : A =(x,y)
 533    if (znC.lt.0.) inC = -inC
        inX = 1
        inY = 2
        goto 534
 534    continue
        cPT1A(1,1) = cPT1C(1,inX)
        cPT1A(1,2) = cPT1C(1,inY)
        cPT2A(1,1) = cPT2C(1,inX)
        cPT2A(1,2) = cPT2C(1,inY)
        if (inC.gt.0) then
          cPT1A(2,1) = cPT1C(2,inX)
          cPT1A(2,2) = cPT1C(2,inY)
          cPT2A(2,1) = cPT2C(2,inX)
          cPT2A(2,2) = cPT2C(2,inY)
          cPT1A(3,1) = cPT1C(3,inX)
          cPT1A(3,2) = cPT1C(3,inY)
          cPT2A(3,1) = cPT2C(3,inX)
          cPT2A(3,2) = cPT2C(3,inY)
        else
          cPT1A(2,1) = cPT1C(3,inX)
          cPT1A(2,2) = cPT1C(3,inY)
          cPT2A(2,1) = cPT2C(3,inX)
          cPT2A(2,2) = cPT2C(3,inY)
          cPT1A(3,1) = cPT1C(2,inX)
          cPT1A(3,2) = cPT1C(2,inY)
          cPT2A(3,1) = cPT2C(2,inX)
          cPT2A(3,2) = cPT2C(2,inY)
        endif
        cPT1A(4,1) = cPT1A(1,1)
        cPT1A(4,2) = cPT1A(1,2)
* Calcul des cotes des triangles T1A et T2A (12,23,31)
* Utile pour connaitre ensuite la normale aux cotes des triangles
        vT1A(1,1) = cPT1A(2,1) - cPT1A(1,1)
        vT1A(1,2) = cPT1A(2,2) - cPT1A(1,2)
        vT1A(2,1) = cPT1A(3,1) - cPT1A(2,1)
        vT1A(2,2) = cPT1A(3,2) - cPT1A(2,2)
        vT1A(3,1) = cPT1A(1,1) - cPT1A(3,1)
        vT1A(3,2) = cPT1A(1,2) - cPT1A(3,2)
C*nu        vT2A(1,1) = cPT2A(2,1) - cPT2A(1,1)
        vT2A(1,2) = cPT2A(2,2) - cPT2A(1,2)
C*nu        vT2A(2,1) = cPT2A(3,1) - cPT2A(2,1)
        vT2A(2,2) = cPT2A(3,2) - cPT2A(2,2)
C*nu        vT2A(3,1) = cPT2A(1,1) - cPT2A(3,1)
        vT2A(3,2) = cPT2A(1,2) - cPT2A(3,2)
* Calcul de la surface de chacun des triangles T1A et T2A
* (en fait, on calcule le double, mais par la suite on ne considere que
*  des rapports de surfaces de triangle)
        SuT1A =  (cPT1A(2,1)+cPT1A(1,1)) * vT1A(1,2)
     &         + (cPT1A(3,1)+cPT1A(2,1)) * vT1A(2,2)
     &         + (cPT1A(1,1)+cPT1A(3,1)) * vT1A(3,2)
        SuT2A =  (cPT2A(2,1)+cPT2A(1,1)) * vT2A(1,2)
     &         + (cPT2A(3,1)+cPT2A(2,1)) * vT2A(2,2)
     &         + (cPT2A(1,1)+cPT2A(3,1)) * vT2A(3,2)
C*DBG        write(ioimp,*) 'Surfaces T1A - T2A :', 0.5*SuT1A,0.5*SuT2A
C*DBG        if (SuT1A.gt.100.*SuT2A .or. SuT2A.gt.100.*SuT1A)
C*DBG     &    write(ioimp,*) 'Rapport des surfaces tres important !'
 
* On initialise l'intersection des triangles avec T2A
        nPIn0 = 3
        do 540 i = 1, nPIn0
          cPIn0(i,1) = cPT2A(i,1)
          cPIn0(i,2) = cPT2A(i,2)
 540    continue
* Determination de l'intersection des triangles T1A et T2A
* en regardant progressivement l'intersection avec les cotes de T1A
* Note : On utilise le fait que les polygones traites sont convexes !
        DO 550 i = 1, 3
          vN1x = -vT1A(i,2)
          vN1y =  vT1A(i,1)
          scal = vN1x * cPT1A(i,1) + vN1y * cPT1A(i,2)
*
          ipos = 0
          ineg = 0
          jpos = 0
          jneg = 0
          do 551 j = 1, nPIn0
            test(j) = (vN1x * cPIn0(j,1) + vN1y * cPIn0(j,2)) - scal
            if (test(j).gt.1.d-10) then
*            if (test(j).gt.0.) then
              ipos = ipos + 1
              if (jneg.ne.0 .and. jpos.eq.0) jpos = j
            else if (test(j).lt.-1.d-10) then
*            else if (test(j).lt.0.) then
              ineg = ineg + 1
              if (jneg.eq.0) jneg = j
            else
              test(j) = 0.
            endif
 551      continue
*
* 1) Cas ou ipos = 0 :
* 1.1) il n'y pas d'intersection (ineg=nPIn0)
* 1.2) l'intersection se limite a 1 point (ineg=nPIn0-1)
* 1.3) l'intersection a 1 segment sur une frontiere (ineg=nPIn0-2)
* La surface correspondante est nulle et il n'y a donc pas de contact !
          if (ipos.eq.0) then
C*DBG       write(ioimp,*) 'Intersection a ',nPIn0-ineg,' points'
            goto 520
          endif
* 2) Cas ou ipos > 0 : Il y a une intersection a calculer
* 2.1) ineg = 0 : Tous les points sont "du bon cote" et l'intersection
*      correspond au polygone de depart !
          if (ineg.eq.0) then
C*DBG            write(ioimp,*) 'Intersection du bon cote - rien a faire'
            goto 550
          endif
* 2.2) ineg > 0 : Il y a deux intersections a determiner car il y a
*      au moins un point du "mauvais cote"
          jpos = max(1,jpos)
          jpr = jpos - 1
          if (jpos.eq.1) jpr = nPIn0
          nPIn = 1
          if (test(jpr).lt.0.) then
            r_z = test(jpos) / (test(jpos) - test(jpr))
            cPIn(nPIn,1) =  cPIn0(jpos,1)
     &                    + r_z*(cPIn0(jpr,1) - cPIn0(jpos,1))
            cPIn(nPIn,2) =  cPIn0(jpos,2)
     &                    + r_z*(cPIn0(jpr,2) - cPIn0(jpos,2))
          else
            cPIn(nPIn,1) = cPIn0(jpr,1)
            cPIn(nPIn,2) = cPIn0(jpr,2)
          endif
          do 552 j = 1, ipos
            nPIn = nPIn + 1
            cPIn(nPIn,1) = cPIn0(jpos,1)
            cPIn(nPIn,2) = cPIn0(jpos,2)
            jpr  = jpos
            jpos = jpos + 1
            if (jpos.gt.nPIn0) jpos = 1
 552      continue
          nPIn = nPIn + 1
          if (test(jpos).lt.0.) then
            r_z = test(jpr) / (test(jpr) - test(jpos))
            cPIn(nPIn,1) =  cPIn0(jpr,1)
     &                    + r_z*(cPIn0(jpos,1) - cPIn0(jpr,1))
            cPIn(nPIn,2) =  cPIn0(jpr,2)
     &                    + r_z*(cPIn0(jpos,2) - cPIn0(jpr,2))
          else
            cPIn(nPIn,1) = cPIn0(jpos,1)
            cPIn(nPIn,2) = cPIn0(jpos,2)
          endif
* Mise a jour cPIn0 pour la suite du traitement des intersections...
          nPIn0 = nPIn
          do 554 j = 1, nPIn0
            cPIn0(j,1) = cPIn(j,1)
            cPIn0(j,2) = cPIn(j,2)
 554      continue
C*DBG-F          write(ioimp,*) 'Intersection ',i,' a ',nPIn0,' points'
C*DBG-F          do j = 1, nPIn0
C*DBG-F            write(ioimp,*) '  ',j,nbpts0+6+j,cPIn0(j,1),cPIn0(j,2)
C*DBG-F          enddo
*
 550    CONTINUE
*
* Calcul de la surface de l'intersection  et de son centre de gravite
        r_z  =  cPIn0(nPIn0,1)*cPIn0(1,2) - cPIn0(1,1)*cPIn0(nPIn0,2)
        SuPIn = r_z
        xGIn = (cPIn0(nPIn0,1)+cPIn0(1,1)) * r_z
        yGIn = (cPIn0(nPIn0,2)+cPIn0(1,2)) * r_z
        do 560 i = 1, nPIn0-1
          j = i + 1
          r_z  =  cPIn0(i,1)*cPIn0(j,2) - cPIn0(j,1)*cPIn0(i,2)
          SuPIn = SuPIn + r_z
          xGIn  = xGIn + (cPIn0(i,1)+cPIn0(j,1)) * r_z
          yGIn  = yGIn + (cPIn0(i,2)+cPIn0(j,2)) * r_z
 560    continue
        if (SuPIn .lt. 1.E-7*max(SuT1A,SuT2A) ) then
C*DBG          write(ioimp,*) 'Intersection a une surface negligeable !'
          goto 520
        endif
        r_z = X1s3 / SupIn
        SuPIn = 0.5 * SupIn
        xGIn = xGIn * r_z
        yGIn = yGIn * r_z
C*DBG        write(ioimp,*) 'Surface Intersection :',SuPIn
*
* Calcul des coordonnees barycentriques du centre de gravite
* pour le triangle T1A
        xpip1 = xGIn - cPT1A(1,1)
        ypip1 = yGIn - cPT1A(1,2)
        xpip2 = xGIn - cPT1A(2,1)
        ypip2 = yGIn - cPT1A(2,2)
        xpip3 = xGIn - cPT1A(3,1)
        ypip3 = yGIn - cPT1A(3,2)
        b1T1 = xpip2 * ypip3 - ypip2 * xpip3
        b2T1 = xpip3 * ypip1 - ypip3 * xpip1
        b3T1 = xpip1 * ypip2 - ypip1 * xpip2
        bt = b1T1 + b2T1 + b3T1
        if (abs(bt-SuT1A) .gt. 1.E-3)
     &    write(ioimp,*) 'Prob bt-SuT1A',bt,SuT1A
        b1T1 = b1T1 / SuT1A
        b2T1 = b2T1 / SuT1A
        b3T1 = b3T1 / SuT1A
*
* pour le triangle T2A
        xpip1 = xGIn - cPT2A(1,1)
        ypip1 = yGIn - cPT2A(1,2)
        xpip2 = xGIn - cPT2A(2,1)
        ypip2 = yGIn - cPT2A(2,2)
        xpip3 = xGIn - cPT2A(3,1)
        ypip3 = yGIn - cPT2A(3,2)
        b1T2 = xpip2 * ypip3 - ypip2 * xpip3
        b2T2 = xpip3 * ypip1 - ypip3 * xpip1
        b3T2 = xpip1 * ypip2 - ypip1 * xpip2
        bt = b1T2 + b2T2 + b3T2
        if (abs(bt-SuT2A) .gt. 1.E-3)
     &    write(ioimp,*) 'Prob bt-SuT2A',bt,SuT2A
        b1T2 = b1T2 / SuT2A
        b2T2 = b2T2 / SuT2A
        b3T2 = b3T2 / SuT2A
*
* Calcul du jeu
        xjeu =   (b1T2*dq1C + b2T2*dq2C + b3T2*dq3C)
     &         - (b1T1*dp1C + b2T1*dp2C + b3T1*dp3C)
C*DBG        write(ioimp,*) 'Jeu =',xjeu
*
*  Ajout d'une relation triangle-triangle
*
        nelri0 = nelri0 + 1
        nelch = nelch + 1
*
*   on ajuste les differents segments si necesaire
*
        if (nelri0.gt.nelrig) then
          nelrig = nelrig + incnel
          segadj,xmatri
          nbelem = nbelem + incnel
          nbnn = 7
          segadj,meleme
          nbnn = 1
          segadj,ipt7
          n = n + incnel
          segadj,mpoval
        endif
*
*       Choix du mult de Lagrange qui porte la condition
*       -> l'element le plus grand
        imult=ipt1.num(1,iel8)
        ielt=iel8
        if (RayT1.lt.RayT2) then
          imult=ipt1.num(1,nbel8+iel6)
          ielt=nbel8+iel6
        endif
*
        num(1,nelri0) = imult
        num(2,nelri0) = ip1
        num(3,nelri0) = ip2
        num(4,nelri0) = ip3
        num(5,nelri0) = jp1
        num(6,nelri0) = jp2
        num(7,nelri0) = jp3
        icolor(nelri0) = 1
*
*  Mise a jour de xmatri
*
*  lambda
        re( 1,1,nelri0) = 0.d0
*  ip1
        re( 2,1,nelri0) =  xnC * b1T1 * SuPIn
        re( 3,1,nelri0) =  ynC * b1T1 * SuPIn
        re( 4,1,nelri0) =  znC * b1T1 * SuPIn
*  ip2
        re( 5,1,nelri0) =  xnC * b2T1 * SuPIn
        re( 6,1,nelri0) =  ynC * b2T1 * SuPin
        re( 7,1,nelri0) =  znC * b2T1 * SuPIn
*  ip3
        re( 8,1,nelri0) =  xnC * b3T1 * SuPIn
        re( 9,1,nelri0) =  ynC * b3T1 * SuPIn
        re(10,1,nelri0) =  znC * b3T1 * SuPIn
*  jp1
        re(11,1,nelri0) = -xnC * b1T2 * SuPIn
        re(12,1,nelri0) = -ynC * b1T2 * SuPIn
        re(13,1,nelri0) = -znC * b1T2 * SuPIn
*  jp2
        re(14,1,nelri0) = -xnC * b2T2 * SuPIn
        re(15,1,nelri0) = -ynC * b2T2 * SuPIn
        re(16,1,nelri0) = -znC * b2T2 * SuPIn
*  jp3
        re(17,1,nelri0) = -xnC * b3T2 * SuPIn
        re(18,1,nelri0) = -ynC * b3T2 * SuPIn
        re(19,1,nelri0) = -znC * b3T2 * SuPIn
*  on transpose
        do 580 ic = 2, nligrp
          re(1,ic,nelri0) = re(ic,1,nelri0)
 580    continue
*  le reste est nul
*
*  remplissage du champoint de depi (jeu)
*
        ipt7.num(1,nelch) = imult
        vpocha(nelch,1)   = (xjeu - xjr) * SuPIn
        vpocha(nelch,2)   = SuPIn
        IF (MELVA2.NE.0) THEN
          NEL2 = min (ielt,NELA)
          VPOCHA(nelch,3) = max(MELVA2.VELCHE(NPTELA,NEL2),0.D0)*SuPIn
        ENDIF
C*DBG-F        call prfaible(ipt1,iel,nPIn0,mfaible)
*
 520  CONTINUE
 519  CONTINUE
*
* Ajustement au plus juste puis desactivation des segments lies
* a la rigidite du type 0
      if (nelri0.ne.nelrig) then
        nelrig = nelri0
        segadj,xmatri
        nbelem = nelri0
        nbnn = 7
        segadj,meleme
      endif
      segdes,xmatri
*
**     write(ioimp,*) ' nb relation type 0 ',nelri0
*      if (nelri0.eq.0) irigel(6,nrigel)=0
*
* Destruction des segments locaux devenus inutiles
      segsup,mfaible
*
C*    GOTO 1000
 
*=======================================================================
*= Fin du sous-programme :
*=======================================================================
 1000 CONTINUE
*
* Ajustement au plus juste du chpoint de depi (jeu) : mpoval et ipt7
* puis desactivation du chpoint
      if (vpocha(/1).ne.nelch) then
        n = nelch
        nc=vpocha(/2)
        segadj,mpoval
        nbnn   = 1
        nbelem = nelch
        nbsous = 0
        nbref  = 0
        segadj,ipt7
**      write(6,*) ' ipt7 dans impo32 '
      endif
* Desctivation de la matrice de raideur de contact
      segdes,mrigid
*
* Reunion des relations portant sur le meme multiplicateur de lagrange
*
      call impofu(MRIGID,MCHPOI)
*     write(6,*) ' apres impofu dans impo32 '
*     call ecchpo(mchpoi,1)
*     call prrigi(mrigid,1)
 
 
*
* Nettoyage eventuel des termes petits dans les relations
*
* A voir plus tard. Frig3c suppose que la relation est complete.
*
**    call relasi(mrigid)
*
 900  CONTINUE
      end