C IMPOS7    SOURCE    PV090527  26/04/30    21:15:41     12529          
      subroutine impos7
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
-INC CCREEL

-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
-INC SMELEME
-INC SMRIGID
-INC SMCHPOI
      dimension vec(3)
      character*4 modepl(4),moforc(4),mnez(3),mcle(3)
      data modepl /'UX  ','UY  ','UR  ','UZ  '/
      data moforc /'FX  ','FY  ','FR  ','FZ  '/
      data mnez   /'PLAT','HEMI','CONE'/
      data mcle   /'MAIT','VECT','ANGL'/
*
*
*
*     lecture des options
*
      call lirmot(mcle,3,ire,0)
*
*     creation de l'objet MAILLAGE contenant tous les contacts possibles
*
      if (ire.eq.1) then
         call impos8
         return
      endif
*
*     lecture du maillage support des conditions
*
      call lirobj('MAILLAGE',ipt1,1,iretou)
      if (ierr.ne.0) return
      call actobj('MAILLAGE',ipt1,1)
      if (ierr.ne.0) return
      if (ipt1.itypel.ne.22) call erreur(851)
      
      SEGACT,MCOORD
*
*     cas du contact ponctuel
*
      if(ipt1.num(/1).eq.3) then
*
*       on test que le dernier noeud est bien identique
*
         jp=ipt1.num(3,1)
         do 80 iel=2,ipt1.num(/2)
            if(ipt1.num(3,iel).ne.jp) call erreur(851)
 80      continue
         call lirmot(mnez,3,iret,0)
*
*       projectile sans nez
*
         if (iret.eq.0) then
            call lirobj('POINT',noeud,1,iretou)
            do 1 ia=1,idim
               vec(ia)=xcoor((noeud-1)*(idim+1)+ia)
 1          CONTINUE
         endif
*
*       projectile a nez plat
*
         if (iret.eq.1) then
            call lirree(vlarg,0,iretou)
            call lirmot(mcle,3,ire,0)
            if(ire.eq.2) then
               call lirobj('POINT',noeud,0,iretou)
               do 2 ia=1,idim
                  vec(ia)=xcoor((noeud-1)*(idim+1)+ia)
 2             CONTINUE
            else
               call erreur(20)
            endif
         endif
*
*       projectile a nez hemispherique
*
         if (iret.eq.2) then
            call lirree(vrayo,0,iretou)
            call lirmot(mcle,3,ire,0)
            if(ire.eq.2) then
               call lirobj('POINT',noeud,0,iretou)
               do 3 ia=1,idim
                  vec(ia)=xcoor((noeud-1)*(idim+1)+ia)
 3             CONTINUE
            else
               call erreur(20)
            endif
         endif
*
*       projectile a nez conique
*
         if (iret.eq.3) then
            call lirree(vlarg,0,iretou)
            call lirmot(mcle,3,ire,0)
            if(ire.eq.3) then
               call lirree(vangl,0,iretou)
               vangl=(XPI/2.D0)-(XPI*vangl/180)
            else
               call erreur(15)
            endif
            call lirmot(mcle,3,ire,0)
            if(ire.eq.2) then
               call lirobj('POINT',noeud,0,iretou)
               do 4 ia=1,idim
                  vec(ia)=xcoor((noeud-1)*(idim+1)+ia)
 4             CONTINUE
            else
               call erreur(20)
            endif
         endif
      endif
*
*     creation de la raideur de conditions
*
      nrige=7
      nrigel=1
      segini mrigid
      ichole=0
      imgeo1=0
      imgeo2=0
      isupeq=0
      iforig=ifour
      coerig(1)=1.d0
      mtymat='RIGIDITE'
*
*     creation du meleme associe a la relation
*
      nbsous=0
      nbref=0
      nbnn=4
      nbelem=0
      segini meleme
      itypel=22
      irigel(1,1)=meleme
*
*     impact ligne-ligne (vrai) ou ligne-noeud (sinon) ?
*
      if(ipt1.num(/1).eq.4) then
*
*       creation des descriptifs communs a toutes les raideurs
*
         imo=1
         if(ifour.eq.0) imo=3
         nligrp=7
         nligrd=nligrp
         segini descr
         irigel(3,1)=descr
         lisinc(1)='LX'
         lisdua(1)='FLX'
         noelep(1)=1
         noeled(1)=1
         do 5 i=2,7,2
            lisinc(i)=modepl(imo)
            lisinc(i+1)=modepl(imo+1)
            lisdua(i)=moforc(imo)
            lisdua(i+1)=moforc(imo+1)
            noelep(i)=(i+2)/2
            noelep(i+1)=noelep(i)
            noeled(i)=noelep(i)
            noeled(i+1)=noelep(i)
 5       continue
*
*       creation du segment imatri
*
         nelrig=50
                 nelri=0
         RIGREL=0
         segini xmatri
         irigel(4,1)=xmatri
*
*       ce qu'on cree est une égalité
*
         irigel(6,1)=0
*
*       ce qu'on cree est symetrique
*
         irigel(7,1)=0
*
*       boucle sur le maillage support
         do 10 iel=1,ipt1.num(/2)
*
*        noeud maitre
*
            ip1=ipt1.num(2,iel)
            xp1=xcoor((ip1-1)*(idim+1)+1)
            yp1=xcoor((ip1-1)*(idim+1)+2)
*
*        deuxieme noeud maitre
*
            ip2=ipt1.num(3,iel)
            xp2=xcoor((ip2-1)*(idim+1)+1)
            yp2=xcoor((ip2-1)*(idim+1)+2)
            xr=xp2-xp1
            yr=yp2-yp1
*
*        noeud esclave
*
            jp=ipt1.num(4,iel)
            xp=xcoor((jp-1)*(idim+1)+1)
            yp=xcoor((jp-1)*(idim+1)+2)
*
*        test de recherche du contact
*
            if(xp2.gt.xp1) then
               if(xp.lt.xp1 .or. xp.gt.xp2) goto 10
            else
               if(xp.lt.xp2 .or. xp.gt.xp1) goto 10
            endif
            if((xp-xp1)*yr-(yp-yp1)*xr.gt.0.d0) goto 10
*
*        determination du noeud defenseur fictif
*
            sr2=xr**2+yr**2
            xa=xp2-xp
            ya=yp2-yp
            alpha=(xa*xr+ya*yr)/sr2
            if(abs(alpha).lt.1D-3) alpha=0.D0
            if(alpha.gt.0.999D0) alpha=1.D0
            if(alpha.lt.-0.999D0) alpha=-1.D0
            sr=sqrt(sr2)
            xr=xr/sr
            yr=yr/sr
*
*        creation du xmatri
*
*            segini xmatri
*            imatri=irigel(4,1)
            nelri=nelri+1
                        if(nelri.gt.nelrig) then
                         nelrig=nelrig+50
             RIGREL=0
             segadj xmatri
                        endif
*            imattt(nelrig)=xmatri
*
*        multiplicateur
*
            re(1,1,nelri)= 0.d0
*
*        1er noeud maitre ip1
*
            re(2,1,nelri)=  yr * alpha
            re(3,1,nelri)= -xr * alpha
*
*        2eme noeud maitre ip2
*
            re(4,1,nelri)=  yr * (1-alpha)
            re(5,1,nelri)= -xr * (1-alpha)
*
*        noeud esclave jp
            re(6,1,nelri)= -yr
            re(7,1,nelri)=  xr
*
*        test de conditionnement
            do 90 il=1,7
               if(abs(re(il,1,nelri)).lt.1D-6) re(il,1,nelri)=0.D0
               if(re(il,1,nelri).gt.(1.D0-1D-6)) re(il,1,nelri)=1.D0
               if(re(il,1,nelri).lt.(-1.D0+1D-6)) re(il,1,nelri)=-1.D0
 90         continue
*
*        transposition du xmatri C
            do 30 il=1,1
               do 301 ic=1,7
                  re(il,ic,nelri)=re(ic,il,nelri)
 301           continue
 30         continue
*
*        modification du meleme
*
            nbelem=num(/2)+1
            nbnn=4
            nbsous=0
            nbref=0
            segadj meleme
            num(1,nbelem)=ipt1.num(1,iel)
            num(2,nbelem)=ip1
            num(3,nbelem)=ip2
            num(4,nbelem)=jp
*
*            segdes xmatri
*
 10      continue
        if(nelri.ne.nelrig) then
                  nelrig=nelri
                  RIGREL=0
                  segadj xmatri
                endif
*
      else
*
*       creation des descriptifs communs a toutes les raideurs
*
         imo=1
         if(ifour.eq.0) imo=3
         nligrp=5
         nligrd=nligrp
         segini descr
         irigel(3,1)=descr
         lisinc(1)='LX'
         lisdua(1)='FLX'
         noelep(1)=1
         noeled(1)=1
         do 70 i=2,5,2
            lisinc(i)=modepl(imo)
            lisinc(i+1)=modepl(imo+1)
            lisdua(i)=moforc(imo)
            lisdua(i+1)=moforc(imo+1)
            noelep(i)=(i+2)/2
            noelep(i+1)=noelep(i)
            noeled(i)=noelep(i)
            noeled(i+1)=noelep(i)
 70      continue
*
*       creation du segment imatri
*
         nelrig=50
                 nelri=0
         RIGREL=0
         segini xmatri
         irigel(4,1)=xmatri
*
*       ce qu'on cree est une égalité
*
         irigel(6,1)=0
*
*       ce qu'on cree est symetrique
*
         irigel(7,1)=0
*
*       boucle sur le maillage support
*
         do 50 iel=1,ipt1.num(/2)
*
*        noeud maitre
*
            ip1=ipt1.num(2,iel)
            xp1=xcoor((ip1-1)*(idim+1)+1)
            yp1=xcoor((ip1-1)*(idim+1)+2)
*
*        deuxieme noeud maitre fictif
*
            if(iel.lt.ipt1.num(/2)) then
               ip2=ipt1.num(2,iel+1)
               jcm=1
            else
               ip2=ipt1.num(2,iel-1)
               jcm=-1
            endif
            xp2=xcoor((ip2-1)*(idim+1)+1)
            yp2=xcoor((ip2-1)*(idim+1)+2)
            xr=(xp2-xp1)*jcm
            yr=(yp2-yp1)*jcm
*
*        noeud esclave
*
            jp=ipt1.num(3,iel)
            xp=xcoor((jp-1)*(idim+1)+1)
            yp=xcoor((jp-1)*(idim+1)+2)
*
*        verification que pas relation du point sur lui meme
*
            if (jp.eq.ip1) goto 50
*
*        detection du contact : point materiel
*
            if (iret.eq.0) then
               if(((xp-xp1)*yr-(yp-yp1)*xr).gt.0.d0) goto 50
               vnorm=sqrt(vec(1)*vec(1)+vec(2)*vec(2))
               vn1=vec(1)/vnorm
               vn2=vec(2)/vnorm
            endif
*
*        detection du contact : projectile a nez plat
*
            if (iret.eq.1) then
               tv=vec(1)/vec(2)
               aa=(xp-xp1)*tv
               if(-(yp+aa-yp1)*xr.gt.0.d0) goto 50
               vnorm=sqrt(vec(1)*vec(1)+vec(2)*vec(2))
               if(abs(xp-xp1).gt.abs(vlarg*vec(2)/vnorm)) goto 50
               vn1=vec(1)/vnorm
               vn2=vec(2)/vnorm
            endif
*
*        detection du contact : projectile a nez hemispherique
*
            if (iret.eq.2) then
               vnorm=sqrt(vec(1)*vec(1)+vec(2)*vec(2))
               if(abs(xp-xp1).gt.vrayo*abs(vec(2))/vnorm) goto 50
               ypo=(vrayo/sqrt(1+(vec(1)*vec(1)/vec(2)*vec(2))))+yp
               xpo=(ypo-yp)*vec(1)/vec(2)+xp
               yy=ypo-sqrt(vrayo*vrayo-(xpo-xp1)*(xpo-xp1))
               if(-(yy-yp1)*xr.gt.0.d0) goto 50
               vnn=sqrt((xpo-xp1)*(xpo-xp1)+(ypo-yp1)*(ypo-yp1))
               vn1=(xp1-xpo)/vnn
               vn2=(yp1-ypo)/vnn
            endif
*
*        detection du contact : projectile a nez conique
*
            if (iret.eq.3) then
               ial=nint(sign(1.D0,(xp-xp1)))
               vangl2=vangl-ial*atan(vec(1)/vec(2))
               aa=abs(xp-xp1)*tan(vangl2)
               if(-(yp+aa-yp1)*xr.gt.0.d0) goto 50
               vnorm=sqrt(vec(1)*vec(1)+vec(2)*vec(2))
               if(abs(xp-xp1).gt.vlarg*abs(vec(2))/vnorm) goto 50
               v1=cos(vangl)*vec(1)+ial*sin(vangl)*vec(2)
               v2=-ial*sin(vangl)*vec(1)+cos(vangl)*vec(2)
               vnn=sqrt(v1*v1+v2*v2)
               vn1=v1/vnn
               vn2=v2/vnn
            endif
*
*        on cree le xmatri C
*
*            segini xmatri
*            imatri=irigel(4,1)
            nelri=nelri+1
                        if(nelri.gt.nelrig) then
                         nelrig=nelrig+50
             RIGREL=0
             segadj xmatri
                        endif
*            imattt(nelrig)=xmatri
*
*        on remplit le xmatri
*
*        multiplicateur
            re(1,1,nelri)= 0.d0
*
*        noeud maitre ip1
            re(2,1,nelri)= vn1
            re(3,1,nelri)=-vn2
*
*        noeud esclave jp
            re(4,1,nelri)=-vn1
            re(5,1,nelri)= vn2
*
*        test de conditionnement
            do 95 il=1,5
               if(abs(re(il,1,nelri)).lt.1D-6) re(il,1,nelri)=0.D0
               if(re(il,1,nelri).gt.(1.D0-1D-6)) re(il,1,nelri)=1.D0
               if(re(il,1,nelri).lt.(-1.D0+1D-6)) re(il,1,nelri)=-1.D0
 95         continue
*
*        transposition du xmatri C
            do 40 il=1,1
               do 401 ic=1,5
                  re(il,ic,nelri)=re(ic,il,nelri)
 401           continue
 40         continue
*
*        modification du meleme
*
            nbelem=num(/2)+1
            nbnn=3
            nbsous=0
            nbref=0
            segadj meleme
            num(1,nbelem)=ipt1.num(1,iel)
            num(2,nbelem)=ip1
            num(3,nbelem)=jp
*
*            segdes xmatri
*
 50      continue
         if(nelrig.ne.nelri) then
                   nelrig=nelri
                   RIGREL=0
                   segadj xmatri
                 endif
*
      endif

      SEGDES,MCOORD
*
*     desactivation des segments
*
*      segdes xmatri,ipt1,mrigid,meleme,descr
*
*     renvoi du resultat
*
      if(nbelem.eq.0) then
         irigel(6,1)=0
         call ecrobj('ANNULE',mrigid)
      else
         call actobj('RIGIDITE',mrigid,1)
         call ecrobj('RIGIDITE',mrigid)
      endif
*
      end