C REDCON    SOURCE    CB215821  24/04/12    21:17:05     11897          
      subroutine redcon(mmodel,mmode2)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)

-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
-INC SMELEME
-INC CCREEL
-INC SMMODEL
      character*4 modepl(3),moforc(3)
      data modepl /'UX  ','UY  ','UZ  '/
      data moforc /'FX  ','FY  ','FZ  '/
*
*  lecture  du maillage support des conditions
*
      segact mmodel
      segini,mmode2=mmodel
      imoin=0
      do  1 imo=1,kmodel(/1)
        imodel=kmodel(imo)
        segact imodel
        do ifo=1,formod(/2)
          if( formod(ifo).eq.'CONTACT' ) go to 2
        enddo
        mmode2.kmodel(imo-imoin)=imodel
        go to 1
    2   continue
* on va reduire le maillage
        ipt1=imamod
        segact ipt1
        segini,ipt2=ipt1
        igard=0
        nbelem=ipt1.num(/2)
        if(idim.eq.3) then
*
*  boucle sur le maillage support en 3D
*

          do 10 iel=1,nbelem
            ip1=ipt1.num(2,iel)
            ip2=ipt1.num(3,iel)
            ip3=ipt1.num(4,iel)
            xp1=xcoor((ip1-1)*(idim+1)+1)
            yp1=xcoor((ip1-1)*(idim+1)+2)
            zp1=xcoor((ip1-1)*(idim+1)+3)
            xp2=xcoor((ip2-1)*(idim+1)+1)
            yp2=xcoor((ip2-1)*(idim+1)+2)
            zp2=xcoor((ip2-1)*(idim+1)+3)
            xp3=xcoor((ip3-1)*(idim+1)+1)
            yp3=xcoor((ip3-1)*(idim+1)+2)
            zp3=xcoor((ip3-1)*(idim+1)+3)
*
*  verification si autre point dans la zone de selection
            jp=ipt1.num(5,iel)
*      write (6,*) 'mail ',ip1,ip2,ip3,jp
            xp =xcoor((jp-1)*(idim+1)+1)
            yp =xcoor((jp-1)*(idim+1)+2)
            zp =xcoor((jp-1)*(idim+1)+3)
*  verification que pas relation du point sur lui meme
            if (jp.eq.ip1) goto 10
            if (jp.eq.ip2) goto 10
            if (jp.eq.ip3) goto 10
*  verif distance au centre de gravite
            xg=(xp1+xp2+xp3)/3d0
            yg=(yp1+yp2+yp3)/3d0
            zg=(zp1+zp2+zp3)/3d0
            d1=sqrt((xg-xp1)**2+(yg-yp1)**2+(zg-zp1)**2)
            d2=sqrt((xg-xp2)**2+(yg-yp2)**2+(zg-zp2)**2)
            d3=sqrt((xg-xp3)**2+(yg-yp3)**2+(zg-zp3)**2)
            dp=sqrt((xg-xp)**2+(yg-yp)**2+(zg-zp)**2)/2d0
            if (dp.gt.d1.and.dp.gt.d2.and.dp.gt.d3) goto 10
*       write (6,*) 'contact test distance ok',dp,d1,d2,d3
*  verif position par rapport aux cotés
*  cote 1 2
            xv=xp2-xp1
            yv=yp2-yp1
            zv=zp2-zp1
            dv=sqrt((xv**2)+(yv**2)+(zv**2))
            if(.not.(dv.lt.1d-10).and.
     >         .not.(dv.gt.1d-10))
     >        write (6,*) 'prob 1 dans impo32'
            if (abs(dv).le.xpetit) write (6,*) 'prob 1 imp32'
            xv=xv/dv
            yv=yv/dv
            zv=zv/dv
            scal=(xp-xp1)*xv+(yp-yp1)*yv+(zp-zp1)*zv
            xm=xp1+scal*xv
            ym=yp1+scal*yv
            zm=zp1+scal*zv
            scal=(xg-xm)*xv+(yg-ym)*yv+(zg-zm)*zv
            xn=(xg-xm)-scal*xv
            yn=(yg-ym)-scal*yv
            zn=(zg-zm)-scal*zv
            dn=sqrt(xn**2+yn**2+zn**2)
            scal=(xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn
            dpm=sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2)
            if (dpm.gt.1d-3*dv.and.scal .lt.-0.2d0*dn*dpm) goto 10
*       write (6,*) 'contact test position 1 ok',scal/(dn*dpm),scal,
*    >  dn,dpm
*  cote 2 3
            xv=xp3-xp2
            yv=yp3-yp2
            zv=zp3-zp2
            dv=sqrt((xv**2)+(yv**2)+(zv**2))
            if(.not.(dv.lt.1d-10).and.
     >         .not.(dv.gt.1d-10))
     >         write (6,*) 'prob 2 dans impo32'
            if (abs(dv).le.xpetit) write (6,*) 'prob 2 imp32'
            xv=xv/dv
            yv=yv/dv
            zv=zv/dv
            scal=(xp-xp2)*xv+(yp-yp2)*yv+(zp-zp2)*zv
            xm=xp2+scal*xv
            ym=yp2+scal*yv
            zm=zp2+scal*zv
            scal=(xg-xm)*xv+(yg-ym)*yv+(zg-zm)*zv
            xn=(xg-xm)-scal*xv
            yn=(yg-ym)-scal*yv
            zn=(zg-zm)-scal*zv
            dn=sqrt(xn**2+yn**2+zn**2)
            scal=(xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn
            dpm=sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2)
            if (dpm.gt.1d-3*dv.and.scal .lt.-0.2d0*dn*dpm) goto 10
*          write (6,*) 'contact test position 2 ok',scal/(dn*dpm),scal,
*    >  dn,dpm
*  cote 3 1
            xv=xp1-xp3
            yv=yp1-yp3
            zv=zp1-zp3
            dv=sqrt((xv**2)+(yv**2)+(zv**2))
            if(.not.(dv.lt.1d-10).and.
     >         .not.(dv.gt.1d-10))
     >         write (6,*) 'prob 3 dans impo32'
            if (abs(dv).le.xpetit) write (6,*) 'prob 3 imp32'
            xv=xv/dv
            yv=yv/dv
            zv=zv/dv
            scal=(xp-xp3)*xv+(yp-yp3)*yv+(zp-zp3)*zv
            xm=xp3+scal*xv
            ym=yp3+scal*yv
            zm=zp3+scal*zv
            scal=(xg-xm)*xv+(yg-ym)*yv+(zg-zm)*zv
            xn=(xg-xm)-scal*xv
            yn=(yg-ym)-scal*yv
            zn=(zg-zm)-scal*zv
            dn=sqrt(xn**2+yn**2+zn**2)
            scal=(xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn
            dpm=sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2)
            if (dpm.gt.1d-3*dv.and.scal .lt.-0.2d0*dn*dpm) goto 10
*       write (6,*) 'contact test position 3 ok',scal/(dn*dpm),scal,
*    >  dn,dpm
* on a un point ou imposer la relation
* Quelle est la relation ???
*     write (6,*) 'contact potentiel '
*
*   direction de la relation
*

            xn=(yp1-yp2)*(zp2-zp3)-(zp1-zp2)*(yp2-yp3)
            yn=(zp1-zp2)*(xp2-xp3)-(xp1-xp2)*(zp2-zp3)
            zn=(xp1-xp2)*(yp2-yp3)-(yp1-yp2)*(xp2-xp3)
            dn=sqrt(xn**2+yn**2+zn**2)
            if(.not.(dn.lt.1d-10).and.
     >         .not.(dn.gt.1d-10))
     >         write (6,*) 'prob 4 dans impo32'
            if (abs(dn).le.xpetit) write (6,*) 'prob 4 imp32'
            xn=xn/dn
            yn=yn/dn
            zn=zn/dn
            if (abs(xn).lt.1d-10) xn=0.d0
            if (abs(yn).lt.1d-10) yn=0.d0
            if (abs(zn).lt.1d-10) zn=0.d0
*
*   recherche du point a mettre en relation (coordonnees barycentriques)
*
            scal=(xp-xg)*xn+(yp-yg)*yn+(zp-zg)*zn
*
*  verif bon cote (a la tolerance pres)
*
            if (scal.lt.-0.2d0*d1.and.scal.lt.-0.2d0*d2.and.
     >          scal.lt.-0.2d0*d3) goto 10

*
*   on range dans le meleme
*

            igard=igard+1
            ipt2.num(1,igard)=ipt1.num(1,iel)
            ipt2.num(2,igard)=ip1
            ipt2.num(3,igard)=ip2
            ipt2.num(4,igard)=ip3
            ipt2.num(5,igard)=jp
   10     continue
          segdes imodel,ipt1
          if(igard.eq.nbelem) then
            mmode2.kmodel(imo-imoin)=imodel
            segsup ipt2
            go to 1
          endif
          if( igard.eq.0) then
            imoin=imoin+1
            segsup ipt2
            go to 1
          endif
          nbelem=igard
          nbnn=ipt2.num(/1)
          nbref=0
          nbsous=0
          segadj ipt2
          segini,imode2=imodel
          imode2.imamod=ipt2
          kmodel(imo-imoin)=imode2
          segdes imode2,ipt2
          go to 1
        elseif(idim.eq.2) then
          if( ipt1.num(/1).eq.5) then
*  formulation faible
            do 210 iel=1,ipt1.num(/2)
              ip1=ipt1.num(2,iel)
              ip2=ipt1.num(3,iel)
              ip3=ipt1.num(4,iel)
              ip4=ipt1.num(5,iel)
              xp1=xcoor((ip1-1)*(idim+1)+1)
              yp1=xcoor((ip1-1)*(idim+1)+2)
              xp2=xcoor((ip2-1)*(idim+1)+1)
              yp2=xcoor((ip2-1)*(idim+1)+2)
              xp3=xcoor((ip3-1)*(idim+1)+1)
              yp3=xcoor((ip3-1)*(idim+1)+2)
              xp4=xcoor((ip4-1)*(idim+1)+1)
              yp4=xcoor((ip4-1)*(idim+1)+2)
*
*   critere d'acceptation de l'éléments:
*   une extrémité d'un segment dans le cercle de sélection de l'autre segment
*   orientations respectives correctes des 2 segments
*
*  verification (provisoire) que pas de noeuds doubles.
              if (ip1.eq.ip3) goto 210
              if (ip1.eq.ip4) goto 210
              if (ip2.eq.ip3) goto 210
              if (ip2.eq.ip4) goto 210
*
              xl13=sqrt((xp3-xp1)**2+(yp3-yp1)**2)
              xl14=sqrt((xp4-xp1)**2+(yp4-yp1)**2)
              xl23=sqrt((xp3-xp2)**2+(yp3-yp2)**2)
              xl24=sqrt((xp4-xp2)**2+(yp4-yp2)**2)
              sca312=(xp3-xp1)*(xp3-xp2)+(yp3-yp1)*(yp3-yp2)
              sca412=(xp4-xp1)*(xp4-xp2)+(yp4-yp1)*(yp4-yp2)
              sca134=(xp1-xp3)*(xp1-xp4)+(yp1-yp3)*(yp1-yp4)
              sca234=(xp2-xp3)*(xp2-xp4)+(yp2-yp3)*(yp2-yp4)
              if (sca312/(xl13*xl23).gt.0.7.and.
     >            sca412/(xl14*xl24).gt.0.7.and.
     >            sca134/(xl13*xl14).gt.0.7.and.
     >            sca234/(xl23*xl24).gt.0.7)   goto 210
*
              scal=(xp1-xp2)*(xp3-xp4)+(yp1-yp2)*(yp3-yp4)
              if (scal.gt.0.) goto 210

* on a un element ou imposer la relation
* Quelle est la relation ???
*
*   direction du contact
*
              xr=(xp2-xp1+xp3-xp4)
              yr=(yp2-yp1+yp3-yp4)
              sr2 = xr**2 + yr**2
              sr = sqrt (sr2)
              xr=xr/sr
              yr=yr/sr
*
*   projection des points sur la direction du contact
*
              xl1=xp1*xr+yp1*yr
              xl2=xp2*xr+yp2*yr
              xl3=xp3*xr+yp3*yr
              xl4=xp4*xr+yp4*yr
*
*   distance des points à leurs projections
*
              d1=(xp1-xl1*xr)*yr-(yp1-xl1*yr)*xr
              d2=(xp2-xl2*xr)*yr-(yp2-xl2*yr)*xr
              d3=(xp3-xl3*xr)*yr-(yp3-xl3*yr)*xr
              d4=(xp4-xl4*xr)*yr-(yp4-xl4*yr)*xr
*
*   calcul de l'intersection des projections.
*
              xm1=min(xl1,xl2)
              xm2=max(xl1,xl2)
              xm3=min(xl3,xl4)
              xm4=max(xl3,xl4)
*
*  critere
*
              xcr=max(xm2-xm1,xm4-xm3)*1d-4
              xi=max(xm1,xm3)
              xj=min(xm2,xm4)
*
              if (xi.ge.xj-xcr) goto 210

*   on range dans le meleme
*
              igard=igard+1
              ipt2.num(1,igard)=ipt1.num(1,iel)
              ipt2.num(2,igard)=ip1
              ipt2.num(3,igard)=ip2
              ipt2.num(4,igard)=ip3
              ipt2.num(5,igard)=ip4
*
  210       continue
            segdes imodel,ipt1
            if(igard.eq.nbelem) then
              mmode2.kmodel(imo-imoin)=imodel
              segsup ipt2
              go to 1
            endif
            if( igard.eq.0) then
              imoin=imoin+1
              segsup ipt2
              go to 1
            endif
            nbelem=igard
            nbnn=ipt2.num(/1)
            nbref=0
            nbsous=0
            segadj ipt2
            segini,imode2=imodel
            imode2.imamod=ipt2
            kmodel(imo-imoin)=imode2
            segdes imode2,ipt2
            go to 1
          elseif(ipt1.num(/1).eq.4) then
* formulation 1 point - 1 segment
            do 110 iel=1,ipt1.num(/2)
              ip1=ipt1.num(2,iel)
              ip2=ipt1.num(3,iel)
              xp1=xcoor((ip1-1)*(idim+1)+1)
              yp1=xcoor((ip1-1)*(idim+1)+2)
              xp2=xcoor((ip2-1)*(idim+1)+1)
              yp2=xcoor((ip2-1)*(idim+1)+2)
*
*  verification si autre point dans le cercle de selection
              jp=ipt1.num(4,iel)
*  verification que pas relation du point sur lui meme
              if (jp.eq.ip1) goto 110
              if (jp.eq.ip2) goto 110
              xp=xcoor((jp-1)*(idim+1)+1)
              yp=xcoor((jp-1)*(idim+1)+2)
              scal=(xp-xp1)*(xp-xp2)+(yp-yp1)*(yp-yp2)
              xl1=sqrt((xp-xp1)**2+(yp-yp1)**2)
              xl2=sqrt((xp-xp2)**2+(yp-yp2)**2)
              scal=scal/(xl1*xl2)
              if (scal.gt.0.7) goto 110
              igard=igard+1
              ipt2.num(1,igard)=ipt1.num(1,iel)
              ipt2.num(2,igard)=ip1
              ipt2.num(3,igard)=ip2
              ipt2.num(4,igard)=jp
  110       continue
            segdes imodel,ipt1
            if(igard.eq.nbelem) then
              mmode2.kmodel(imo-imoin)=imodel
              segsup ipt2
              go to 1
            endif
            if( igard.eq.0) then
              imoin=imoin+1
              segsup ipt2
              go to 1
            endif
            nbelem=igard
            nbnn=ipt2.num(/1)
            nbref=0
            nbsous=0
            segadj ipt2
            segini,imode2=imodel
            imode2.imamod=ipt2
            mmode2.kmodel(imo-imoin)=imode2
            segdes imode2,ipt2
            go to 1
          endif
        endif
    1 continue
      segdes mmodel
      if(imoin.ne.0) then
         n1=kmodel(/1)-imoin
         if(n1.eq.0) then
            call erreur (21)
            return
         endif
         segadj mmode2
      endif
      segdes mmode2
      return
      end