C IMPOS1    SOURCE    FANDEUR   12/03/15    21:26:25     7314           

C  Operateur 'IMPO' 'MAIL' ( | 'MESC' | ) Lig1 Lig2 ( Crit) ;
C    2D                      | 'FAIB' |
C                            | 'SYME' |

      SUBROUTINE IMPOS1

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

-INC CCOPTIO

-INC SMCOORD
-INC SMELEME

      PARAMETER ( X1s3 = 0.33333333333333333333333333333333333333333D0 ,
     &            X1s4 = 0.25000000000000000000000000000000000000000D0 )

      logical faible
      PARAMETER (ncle = 3)
      character*4 mcle(ncle)
      data mcle / 'MESC','FAIB','SYME' /
*
      idimp1 = IDIM + 1
*
*  Lecture du mot eventuel SYME(trique) ou MESC ou FAIB(le)
*
      mesc = 0
      call lirmot(mcle,ncle,mesc,0)
      if (ierr.ne.0) return
*
*  Lecture des deux lignes en contact
*
      call lirobj('MAILLAGE',ipt1,1,iretou)
      if (ierr.ne.0) return
      call lirobj('MAILLAGE',ipt2,1,iretou)
      if (ierr.ne.0) return
*
*  Lecture eventuelle de la distance limite (crit)
*
      crit=1.e+35
      call lirree(crit,0,iretou)
      if (ierr.ne.0) return
      cri2 = crit*crit
*      crit = abs(crit)
*DBG     write(ioimp,*) ' critere ',crit,cri2
*
*  Petites verifications sur les lignes
*
      segact,ipt1,ipt2
      if (ipt1.itypel.ne.2) call erreur(16)
      if (ipt2.itypel.ne.2) call erreur(16)
      if (ierr.ne.0) goto 900
*
*  Quelques initialisations
*
      nbpts0 = xcoor(/1) / idimp1
*
      nel1 = ipt1.num(/2)
      nel2 = ipt2.num(/2)
*      noe1 = ipt1.num(/1)
      noe2 = ipt2.num(/1)
* En fait : noe2 = noe1 = 2 car elements de type SEG2 !
*
*  Increment sur le nombre d'elements de contact retenus pour agrandir 
*  les longueurs des segments (maillage, mcoord) de maniere adequate
      incnel = 10000
*
      faible = mesc.eq.2
      if (mesc.ne.1) mesc = 0

*=======================================================================
* Formulation forte ("standard") du contact
*=======================================================================
      if (faible) goto 200
*
*  Creation du maillage de contact (meleme)
*
      nbsous = 0
      nbref  = 0
      nbnn   = 4
      nbelem = incnel
      segini,meleme
      itypel = 22
*
*  Nombre de points supports (multiplicateurs de Lagrange) associes
*  au contact : 1 par element de contact
*  Ajustement du segment mcoord
*
      nbpts = nbpts0 + (incnel * 1)
      segadj,mcoord
*
      nblag = nbpts0
      nbcoo = (nblag-1)*idimp1
      nbelc = 0
*
      ipas = 1
*
 5    continue
*
*  On boucle sur les elements de la 1ere ligne de contact
*
      do 10 iel = 1, nel1
        nbeini = nbelc
        ip1=ipt1.num(1,iel)
        ip2=ipt1.num(2,iel)
        xp1=xcoor((ip1-1)*idimp1+1)
        yp1=xcoor((ip1-1)*idimp1+2)
        xp2=xcoor((ip2-1)*idimp1+1)
        yp2=xcoor((ip2-1)*idimp1+2)
        xg1 = xp1+xp2
        yg1 = yp1+yp2
*
*  Boucle sur les noeuds de la 2eme ligne de contact
*
        do 20 jel=1,nel2
          do 20 jnn=1,noe2
            jp=ipt2.num(jnn,jel)
*  verification que pas relation du point sur lui meme
            if (jp.eq.ip1) goto 20
            if (jp.eq.ip2) goto 20
*  verification que pas deja la relation
            do 21 irela = nbeini+1, nbelc
              if (jp.eq.num(4,irela)) goto 20
 21         continue
            xp=xcoor((jp-1)*idimp1+1)
            yp=xcoor((jp-1)*idimp1+2)
*  test distance
            d1 = (xp1-xp)**2+(yp1-yp)**2
            d2 = (xp2-xp)**2+(yp2-yp)**2
            if (cri2.lt.d1.and.cri2.lt.d2) goto 20
*
*  Mise a jour de la dimension des segments si necessaire
*
            if (nbelc.eq.nbelem) then
              nbelem = nbelem + incnel
              segadj,meleme
              nbpts  = nbpts  + incnel
              segadj,mcoord
            endif
*
*   on incremente mcoord avec le pt support du multiplicateur
*   on range dans le meleme l'element de contact associe
*
            nblag = nblag + 1
            nbcoo = nbcoo + idimp1
            nbelc = nbelc + 1
*
            xcoor(nbcoo+1) = (xg1+xp) * X1s3
            xcoor(nbcoo+2) = (yg1+yp) * X1s3
            xcoor(nbcoo+3) = 0.
*
            num(1,nbelc) = nblag
            num(2,nbelc) = ip1
            num(3,nbelc) = ip2
            num(4,nbelc) = jp
*
20      continue
*
10    continue
*
*  On recommence en inversant les deux lignes en cas de traitement 
*  symetrique du contact (option MESC non activee)
*
      if (ipas.eq.1.and.mesc.ne.1) then
        ipas=ipas+1
        iaux=ipt1
        ipt1=ipt2
        ipt2=iaux
        nel1 = ipt1.num(/2)
        nel2 = ipt2.num(/2)
        noe2 = ipt2.num(/1)
        goto 5
      endif
*
      goto 1000

*=======================================================================
* Formulation "faible" du contact
*=======================================================================
 200  continue
*
*  Le nb d'elements maximal a creer pour le maillage de contact FAIBLE :
*  1 entre chaque element de chacune des lignes
*
*  Creation du maillage de contact (meleme)
*
      nbsous= 0
      nbref = 0
      nbnn  = 5
      nbelem = incnel
      segini,meleme
      itypel = 22
*
*  Nomre de points supports (multiplicateurs de Lagrange) associes
*  au contact : 1 par element de contact
*  Ajustement du segment mcoord
*
      nbpts = nbpts0 + (nbelem * 1)
      segadj,mcoord
*
      nblag = nbpts0
      nbcoo = (nblag-1)*idimp1
      nbelc = 0
*
*  Boucle sur les elements de la 1ere ligne de contact
*
      do 210 iel = 1, nel1
        ip1 = ipt1.num(1,iel)
        ip2 = ipt1.num(2,iel)
        ipv = (ip1-1)*idimp1
        xp1 = xcoor(ipv+1)
        yp1 = xcoor(ipv+2)
        ipv = (ip2-1)*idimp1
        xp2 = xcoor(ipv+1)
        yp2 = xcoor(ipv+2)
        xg1 = xp1 + xp2
        yg1 = yp1 + yp2
*
*  Boucle sur les elements du maillage de la 2e ligne de contact
*
        do 220 jel=1,nel2
          jp1 = ipt2.num(1,jel)
          jp2 = ipt2.num(2,jel)
          ipv = (jp1-1)*idimp1
          xp3 = xcoor(ipv+1)
          yp3 = xcoor(ipv+2)
          ipv = (jp2-1)*idimp1
          xp4 = xcoor(ipv+1)
          yp4 = xcoor(ipv+2)

*  test distance
          d1 = (xp1-xp3)**2+(yp1-yp3)**2
          d2 = (xp2-xp3)**2+(yp2-yp3)**2
          if (d1.le.cri2.or.d2.le.cri2) goto 250
          d1 = (xp1-xp4)**2+(yp1-yp4)**2
          d2 = (xp2-xp4)**2+(yp2-yp4)**2
          if (d1.le.cri2.or.d2.le.cri2) goto 250
          goto 220
*
 250      continue
*
*  Mise a jour de la dimension des segments si necessaire
*
          if (nbelc.eq.nbelem) then
            nbelem = nbelem + incnel
            segadj,meleme
            nbpts  = nbpts  + incnel
            segadj,mcoord
          endif
*
*   on incremente mcoord avec les pts support des multiplicateurs
*   on range dans le meleme l'element de contact correspondant
*
          nblag = nblag + 1
          nbcoo = nbcoo + idimp1
          nbelc = nbelc + 1
*
          xcoor(nbcoo+1) = (xg1+xp3+xp4) * X1s4
          xcoor(nbcoo+2) = (yg1+yp3+yp4) * X1s4
          xcoor(nbcoo+3) = 0.
*
          num(1,nbelc) = nblag
          num(2,nbelc) = ip1
          num(3,nbelc) = ip2
          num(4,nbelc) = jp1
          num(5,nbelc) = jp2
*
 220    continue
*
 210  continue
*
*     goto 1000
*
*=======================================================================
* Fin
*=======================================================================
 1000 continue
*DBG      write(ioimp,*) ' nbpts nbpts0 nblag',nbpts,nbpts0,nblag
*DBG      write(ioimp,*) ' nbelem nbelc',nbelem,nbelc
*
*  On renvoie le resultat mais au prealable on ajuste les longueurs
*
      if (nbelem.lt.nbelc) then
        call erreur(5)
        segsup,meleme
        goto 900
      else if (nbelem.gt.nbelc) then
        nbelem = nbelc
        segadj,meleme
        nbpts = nblag
        segadj,mcoord
      endif
      segdes,meleme
      call ecrobj('MAILLAGE',meleme)

 900  continue
      segdes,ipt1,ipt2

      return
      end