C RELAMI    SOURCE    PV090527  26/04/30    21:16:09     12529          

      SUBROUTINE RELAMI(IP0,IPT1,IPRIG)

C=======================================================================
C CE SOUS-PROGRAMME CONSTRUIT LA RIGIDITE LIANT LINEAIREMENT LES DDL
C DES NOEUDS MILIEUX D'UN MAILLAGE QUADRATIQUE AUX NOEUDS SOMMETS
C
C SYNTHAXE GIBIANE : RIG1 = RELA MILI (LMOTS1) GEO1 ;
C=======================================================================

C======================== ZONE DE DECLARATIONS =========================

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

-INC PPARAM
-INC CCOPTIO
-INC CCHAMP
-INC CCGEOME
-INC SMLMOTS
-INC SMCOORD
-INC SMELEME
-INC SMRIGID
-INC CCREEL 

      CHARACTER*4 LESDDL(10),LESDUA(10)

      SEGMENT IMILI(NBDDL)

C========================= CORPS DU PROGRAMME ==========================

C==== NOMS DES INCONNUES DE LA RIGIDITE

C     Si pas de LISTMOTS, je prends le inconnues de la meca.
      IF (IP0.EQ.0) THEN
        IF (IFOUR.LT.0) THEN
          NBDDL=2
          LESDDL(1)='UX  '
          LESDDL(2)='UY  '
          LESDUA(1)='FX  '
          LESDUA(2)='FY  '
        ELSEIF (IFOUR.EQ.0) THEN
          NBDDL=2
          LESDDL(1)='UR  '
          LESDDL(2)='UZ  '
          LESDUA(1)='FR  '
          LESDUA(2)='FZ  '
        ELSEIF (IFOUR.EQ.1) THEN
          NBDDL=3
          LESDDL(1)='UR  '
          LESDDL(2)='UT  '
          LESDDL(3)='UZ  '
          LESDUA(1)='FR  '
          LESDUA(2)='FT  '
          LESDUA(3)='FZ  '
        ELSE
          NBDDL=3
          LESDDL(1)='UX  '
          LESDDL(2)='UY  '
          LESDDL(3)='UZ  '
          LESDUA(1)='FX  '
          LESDUA(2)='FY  '
          LESDUA(3)='FZ  '
        ENDIF
      ELSE
C     Sinon, on prends les DDL specifies et on cherche les duals
        MLMOTS=IP0
        SEGACT,MLMOTS
        NBDDL=MOTS(/2)
        DO I=1,NBDDL
           DO J=1,LNOMDD
             IF (MOTS(I).EQ.NOMDD(J)) THEN
               LESDDL(I)=NOMDD(J)
               LESDUA(I)=NOMDU(J)
             ENDIF
           ENDDO
        ENDDO
        SEGDES,MLMOTS
      ENDIF

C==== TRANSFORMATION DU MAILLAGE INI. EN SEGMENTS SI BESOIN

      IF (IDIM.GE.2) THEN
        CALL ECROBJ('MAILLAGE',IPT1)
        CALL CHANLG
        IF (IPT1.EQ.0) THEN
          CALL ERREUR(16)
        END IF
        CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU)
        IF (IERR.NE.0) RETURN
      ENDIF

C==== CONSTRUCTION DU MAILLAGE SUPPORT DE LA RIGIDITE

C     J'initialise un vecteur que je vais remplir de maillages elem.
      SEGINI,IMILI

      IDIMP1=IDIM+1
C     Je parcours le maillage ini. et construis les maillages elem.
      SEGACT,IPT1
      NBSOUS1=IPT1.LISOUS(/1)
C     J'ai un maillage simple
      IF (NBSOUS1.EQ.0) THEN
        IF (IPT1.ITYPEL.EQ.3) THEN
          NBEL1=IPT1.NUM(/2)
          segact mcoord*mod
          NBPTI=nbpts
          NBPTS=NBPTI+NBDDL*NBEL1
          SEGADJ,MCOORD
          NBSOUS=0
          NBREF=0
          NBELEM=NBEL1
          NBNN=4
          ICPT1=0
          DO I=1,NBDDL
            SEGINI,IPT2
            IPT2.ITYPEL=22
            DO J=1,NBEL1
              IP=NBPTI+ICPT1+J
              IREF=(IP-1)*IDIMP1
              IREF2=(IPT1.NUM(2,J)-1)*IDIMP1
              DO K=1,IDIMP1
                XCOOR(IREF+K)=XCOOR(IREF2+K)
              ENDDO
              IPT2.NUM(1,J)=IP
              IPT2.NUM(2,J)=IPT1.NUM(2,J)
              IPT2.NUM(3,J)=IPT1.NUM(1,J)
              IPT2.NUM(4,J)=IPT1.NUM(3,J)
            ENDDO
            SEGDES,IPT2
            IMILI(I)=IPT2
            ICPT1=ICPT1+NBEL1
          ENDDO
        ELSE
C         Si pas de SEG3, ERREUR
          SEGSUP,IMILI
          SEGDES,IPT1
          CALL ERREUR(16)
          RETURN
        ENDIF
C     J'ai un maillage complexe
      ELSE
        NBEL1=0
        DO ISOUS=1,NBSOUS1
          IPT3=IPT1.LISOUS(ISOUS)
          SEGACT,IPT3
          IF (IPT3.ITYPEL.EQ.3) THEN
            NBEL1=NBEL1+IPT3.NUM(/2)
          ENDIF
        ENDDO
C       Si pas de SEG3, ERREUR
        IF (NBEL1.EQ.0) THEN
          SEGSUP,IMILI
          SEGDES,IPT1
          CALL ERREUR(16)
          RETURN
        ENDIF
        segact mcoord*mod
        NBPTI=nbpts
        NBPTS=NBPTI+NBDDL*NBEL1
        SEGADJ,MCOORD
        NBSOUS=0
        NBREF=0
        NBELEM=NBEL1
        NBNN=4
        ICPT1=0
        DO I=1,NBDDL
          ICEL1=0
          SEGINI,IPT2
          IPT2.ITYPEL=22
          DO ISOUS=1,NBSOUS1
            IPT3=IPT1.LISOUS(ISOUS)
            IF (IPT3.ITYPEL.EQ.3) THEN
              NBEL3=IPT3.NUM(/2)
              DO J=1,NBEL3
                IP=NBPTI+ICPT1+J
                IREF=(IP-1)*IDIMP1
                IREF2=(IPT3.NUM(2,J)-1)*IDIMP1
                DO K=1,IDIMP1
                  XCOOR(IREF+K)=XCOOR(IREF2+K)
                ENDDO
                IPT2.NUM(1,ICEL1+J)=IP
                IPT2.NUM(2,ICEL1+J)=IPT3.NUM(2,J)
                IPT2.NUM(3,ICEL1+J)=IPT3.NUM(1,J)
                IPT2.NUM(4,ICEL1+J)=IPT3.NUM(3,J)
              ENDDO
              ICPT1=ICPT1+NBEL3
              ICEL1=ICEL1+NBEL3
            ENDIF
            SEGDES,IPT3
          ENDDO
          SEGDES,IPT2
          IMILI(I)=IPT2
        ENDDO
      ENDIF
      SEGDES,IPT1

C==== CONSTRUCTION DE LA RIGIDITE ASSOCIEE AUX RELA.

      NRIGEL=NBDDL
      NRIGE=8
      SEGINI,RI1
      RI1.MTYMAT='RIGIDITE'
      RI1.IFORIG=IFOUR
      NLIGRP=4
      NLIGRD=NLIGRP
      DO I=1,NRIGEL
C       On a un segment DESCR par DDL
        SEGINI,DES1
        DES1.LISINC(1)='LX'
        DES1.LISDUA(1)='FLX'
        DES1.NOELEP(1)=1
        DES1.NOELED(1)=1
        DO J=2,4
          DES1.LISINC(J)=LESDDL(I)
          DES1.LISDUA(J)=LESDUA(I)
          DES1.NOELEP(J)=J
          DES1.NOELED(J)=J
        ENDDO
        SEGDES,DES1
C       On a un segment XMATRI par DDL
        nelrig=nbel1
        rigrel=0
        SEGINI,XMATR1
        XMATR1.RE(1,2,1)=1.D0
        XMATR1.RE(1,3,1)=-0.5D0
        XMATR1.RE(1,4,1)=-0.5D0
        XMATR1.RE(2,1,1)=XMATR1.RE(1,2,1)
        XMATR1.RE(3,1,1)=XMATR1.RE(1,3,1)
        XMATR1.RE(4,1,1)=XMATR1.RE(1,4,1)
*        SEGDES,XMATR1
C       On a NBEL1 matrice(s) elementaire(s)
        do ioup=2,nelrig
           do io=1,xmatr1.re(/2)
              do iu=1,xmatr1.re(/1)
                 xmatr1.re(iu,io,ioup)=xmatr1.re(iu,io,1)
               enddo
           enddo
        enddo
*  reprendre xmatr1 pour ponderer en fonction des positions reelles des noeuds.
          meleme=imili(i)
         do ioup=1,nelrig
          segact meleme
          ipm=num(2,ioup)
          ip1=num(3,ioup)
          ip2=num(4,ioup)
          xpm=xcoor((ipm-1)*(idim+1)+1)
          ypm=xcoor((ipm-1)*(idim+1)+2)
          if (idim.eq.3) then              
          zpm=xcoor((ipm-1)*(idim+1)+3)
          else
          zpm=0.d0
          endif
          xp1=xcoor((ip1-1)*(idim+1)+1)
          yp1=xcoor((ip1-1)*(idim+1)+2)
          if (idim.eq.3) then              
          zp1=xcoor((ip1-1)*(idim+1)+3)
          else
          zp1=0.d0
          endif
          xp2=xcoor((ip2-1)*(idim+1)+1)
          yp2=xcoor((ip2-1)*(idim+1)+2)
          if (idim.eq.3) then              
          zp2=xcoor((ip2-1)*(idim+1)+3)
          else
          zp2=0.d0
          endif
          xl12=(xp1-xp2)*(xp1-xp2)+(yp1-yp2)*(yp1-yp2)+
     >    (zp1-zp2)*(zp1-zp2)
          xl2=max(xl2,xpetit)
          xlm1=(xp1-xpm)*(xp1-xp2)+(yp1-ypm)*(yp1-yp2)+
     >    (zp1-zpm)*(zp1-zp2)
          xcoef=xlm1/xl12
          xcoef=min(1.d0,max(0.d0,xcoef))
        XMATR1.RE(1,3,ioup)=-1.d0+xcoef
        XMATR1.RE(1,4,ioup)=-xcoef
        XMATR1.RE(3,1,ioup)=XMATR1.RE(1,3,ioup)
        XMATR1.RE(4,1,ioup)=XMATR1.RE(1,4,ioup)
         enddo
         segdes xmatr1

C       On remplit la rigidite
        RI1.COERIG(I)=1.
        RI1.IRIGEL(1,I)=IMILI(I)
        RI1.IRIGEL(2,I)=0
        RI1.IRIGEL(3,I)=DES1
        RI1.IRIGEL(4,I)=xMATR1
        RI1.IRIGEL(5,I)=NIFOUR
        RI1.IRIGEL(6,I)=0
        RI1.IRIGEL(7,I)=0
        RI1.IRIGEL(8,I)=0
      ENDDO
      SEGDES,RI1
      IPRIG=RI1
      SEGSUP IMILI

      END