C RELA1     SOURCE    PV090527  26/04/30    21:16:06     12529          
      SUBROUTINE RELA1
C
C     RELAZIONE LINEARE TRA DDL
C
C     ICONR  NUMERO OGGETTI MELEME DELLA RELAZIONE
C     ITYESR(NBSOUS) LISTA TIPI ELEMENTI CONTENUTI IN MELEME-N
C     MUNESR(NBSOUS)  NUMERO ELEMENTI IN OGNI SOTTOSTRUTTURA DI MELEME-N
C     LNODSR(NNR) LISTA NODI MELEME-N
C     IPOR1(N)  PUNTATORE SU MWREL1
C     IPOR2(N)  PUNTATORE SU MWREL2
C     IPOR3(N)  PUNTATORE SU MWREL3
C     COEFR(N)  COEFFICIENTI RELAZIONE
C     INCREL(N) IDENTIFICATORE NOME DDL
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)

-INC SMELEME
-INC SMCOORD
-INC SMRIGID
-INC TMTRAV
-INC SMCHPOI

-INC PPARAM
-INC CCOPTIO
-INC CCREEL
-INC CCGEOME
-INC CCHAMP

      SEGMENT /MWGGM1/(IPORE1(0))
      SEGMENT /MWGGM2/(IPORE2(0))
      SEGMENT /MWGGM3/(IPORE3(0))
      SEGMENT /MWGGM4/(INCREL(0))
      SEGMENT /MWGGM5/(COEFR(0)*D)
      SEGMENT /MWREL1/(ITYESR(0)),IEL11.MWREL1,
     +     IEL21.MWREL1,iel01.MWREL1
      SEGMENT /MWREL2/(MUNESR(0)),IEL12.MWREL2,
     +     IEL22.MWREL2,iel02.MWREL2
      SEGMENT /MWREL3/(LNODSR(0)),IEL13.MWREL3,
     +     IEL23.MWREL3
      CHARACTER*4 MOTPV(3)
      CHARACTER*4 MOCORI(7)
      CHARACTER*4 MOTPM(2)
      CHARACTER*4 MOTDDL(3)
      CHARACTER*4 MOTBLO(3)
      CHARACTER*4 MODEPL(8)
      CHARACTER*4 MOROTA(5)
      CHARACTER*4 MODUAL(1)
      CHARACTER*8 MILMOT
      DIMENSION XNOR(3)
      DATA MOTPV /'MINI','MAXI','FROT'/
      DATA MOTPM /'+   ','-   '/
      DATA MOCORI/'CORI','ENSE','ACCR','GLIS','BARY','MILI','TUYA'/
C
      DATA MOTBLO/'DEPL','ROTA','DIRE'/
      DATA MODEPL/'UX  ','UY  ','UZ  ','UR  ','UZ  ','UT  ',
     &     'ALFA','BETA'/
      DATA MOROTA/'RX  ','RY  ','RZ  ','RT  ','RS  '/
      DATA MODUAL/'DUAL'/
c
      iel01=0
      iel02=0
C
C     EST-CE UNE RELATION DE CORPS RIGIDE ?
C

      CALL LIRMOT(MOCORI,7,ICORI,0)
*      write(ioimp,*) 'rela1, icori=',icori
      IF (ICORI.NE.0) THEN
c     BARY
         if (icori.eq.5) then
            call relaba
c     MILI
         ELSE IF (ICORI.EQ.6) THEN
            CALL LIROBJ('LISTMOTS',IP0,0,IRETOU)
            IF (IERR.NE.0) RETURN
            CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU)
            IF (IERR.NE.0) RETURN
            CALL RELAMI(IP0,IPT1,IPRIG)
            CALL ECROBJ('RIGIDITE',IPRIG)
c     TUYA
         elseif(icori.eq.7) then
            call reltuy
c     CORI, ENSE, ACCRO, GLIS
         else
            CALL RELA2(MOCORI(ICORI))
         endif
         RETURN
      ENDIF
C
C     EST-CE UNE CONDITION UNILATERALE ?
C
      NILATE=0
      CALL LIRMOT (MOTPV,3,IPO,0)
      IF(IPO.EQ.1) NILATE=-1
      IF(IPO.EQ.2) NILATE=+1
      IF(IPO.EQ.3) NILATE=+2
C
C     INITIALISATIONS
C
      COEX=1.D0
C
C     Deformations planes ou contraintes planes ou defo. plane gene :
      IF (IFOUR.EQ.-1.OR.IFOUR.EQ.-2.OR.IFOUR.EQ.-3) THEN
         LDEPL=2
         IADEPL=0
         LROTA=1
         IAROTA=2
C     Axisymetrique :
      ELSE IF (IFOUR.EQ.0) THEN
         LDEPL=2
         IADEPL=3
         LROTA=1
         IAROTA=3
C     Fourier :
      ELSE IF (IFOUR.EQ.1) THEN
         LDEPL=3
         IADEPL=3
         LROTA=2
         IAROTA=2
C     Tridimensionnel :
      ELSE IF (IFOUR.EQ.2) THEN
         LDEPL=3
         IADEPL=0
         LROTA=3
         IAROTA=0
C     Massif 1D (IDIM=1) :
      ELSE IF (IFOUR.GE.3.AND.IFOUR.LE.15) THEN
         LDEPL=1
         IADEPL=0
         IF (IFOUR.GE.12) IADEPL=3
         LROTA=0
         IAROTA=0
C     Autres cas :
      ELSE
         LDEPL=0
         IADEPL=0
         LROTA=0
         IAROTA=0
      ENDIF
C
      SEGINI MWGGM1
      SEGINI MWGGM2
      SEGINI MWGGM3
      SEGINI MWGGM4
      SEGINI MWGGM5
*
*     LECTURE EVENTUELLE D'UN MODELLE
*

      CALL LIROBJ('MMODEL',IPOMOD,0,IREMOD)
      IF(IREMOD.NE.0) THEN
         CALL RELMOD(IPOMOD)
         RETURN
      ENDIF
*
*     LECTURE EVENTUELLE D'UN CHPOINT
*
      CALL LIROBJ('CHPOINT',IPOCHP,0,IRECHP)
      IF(IRECHP.NE.0) THEN
        CALL ACTOBJ('CHPOINT',IPOCHP,1)
*  lecture eventuelle mot cle DUAL
         IRECHD=0
         CALL LIRMOT (MODUAL,1,IRECHD,0)
         IF(IRECHD.EQ.1) GOTO 500
*
*     ON MET LE CHPOINT SOUS FORME DE TABLEAU DE TRAVAIL
*
         CALL TRACHP(IPOCHP,MTRAV)
         SEGACT MTRAV
         DO 1990 I=1,IBIN(/1)
            DO 19901 J=1,IBIN(/2)
               IF(IBIN(I,J).EQ.0) GO TO 19901
               CALL PLACE(NOMDD,LNOMDD,IPO,INCO(I))
               if (ipo.eq.0) then
                 CALL PLACE(NOMDU,LNOMDU,IPO,INCO(I))
                 ipo = -ipo
               endif
               IF(IPO.EQ.0) THEN
                  moterr(1:4)=inco(i)
                  SEGSUP MTRAV
                  call erreur(108)
                  GO TO 559
               ENDIF
               SEGINI MWREL1
               SEGINI MWREL2
               SEGINI MWREL3
               IPORE1(**)=MWREL1
               IPORE2(**)=MWREL2
               IPORE3(**)=MWREL3
               COEFR(**)=BB(I,J)
               INCREL(**)=IPO
               ITYESR(**)=1
               MUNESR(**)=1
               LNODSR(**)=IGEO(J)
19901       CONTINUE
 1990    CONTINUE
         SEGSUP MTRAV
         GO TO 300
      ENDIF
 200  CONTINUE
C
C     LETTURA OPERANDI
C
      SEGINI MWREL1
      SEGINI MWREL2
      SEGINI MWREL3
      IPORE1(**)=MWREL1
      IPORE2(**)=MWREL2
      IPORE3(**)=MWREL3
 1320 CONTINUE
      CALL LIRREE(XR,0,IRETOF)
      IF(IRETOF.EQ.0) THEN
         CALL LIRENT(IR,0,IRETOI)
         XR=1.D0
         IF(IRETOI.NE.0) XR=DBLE(IR)
      ENDIF
      CALL LIRMOT(NOMDD,LNOMDD,IPO,0)
      IF(IPO.NE.0) THEN
         COEFR(**)=XR*COEX
         INCREL(**)=IPO
      ELSE
*
*     ON REGARDE SI IL Y A DEPL OU ROTA SUIVI DE DIRECTION
*
C     En DIMENSION 1, le mot-cle 'DIRE' est interdit.
         IF (IDIM.EQ.1) THEN
            INTERR(1)=IDIM
            MOTERR(1:4)=MOTBLO(3)
            CALL ERREUR(971)
            GOTO 559
         ENDIF
         CALL LIRMOT(MOTBLO,2,IMOT,1)
         IF (IMOT.EQ.0) GOTO 559
         IF (IMOT.EQ.1) THEN
            IBDDL=LDEPL
            DO 4481 IA=1,IBDDL
               MOTDDL(IA)=MODEPL(IADEPL+IA)
 4481       CONTINUE
         ENDIF
         IF (IMOT.EQ.2) THEN
            IBDDL=LROTA
            DO 4482 IA=1,IBDDL
               MOTDDL(IA)=MOROTA(IAROTA+IA)
 4482       CONTINUE
         ENDIF
         CALL LIRMOT(MOTBLO(3),1,IMOT,1)
         IF(IMOT.EQ.0) GO TO 559
         CALL LIROBJ('POINT',KPOINT,1,IRETOU)
         IF(IRETOU.EQ.0) GO TO 559
         YL=0.D0
         XPREC=SQRT(XPETIT/XZPREC)
         segact mcoord
         DO 4484 IA=1,IDIM
            XVAL    =XCOOR((KPOINT-1)*(IDIM+1)+IA)
            XNOR(IA)=XVAL
            IF(ABS(XVAL) .GT. XPREC)THEN
              YL=YL+XVAL**2
            ENDIF
 4484    CONTINUE
         IF (YL .EQ. 0.D0) THEN
            CALL ERREUR(239)
            GOTO 559
         ENDIF
         YL=1.D0/SQRT(YL)
         DO 4485 IA=1,IDIM
            XNOR(IA)=XNOR(IA)*YL
 4485    CONTINUE
*
*     ON LIT LE MAILLAGE
*
         CALL LIROBJ('POINT   ',MILPOI,0,IRETOU)
         IF(IRETOU.NE.0) THEN
            MILMOT='POINT   '
         ELSE
            CALL LIROBJ('MAILLAGE',MILPOI,1,IRETOU)
            IF(IRETOU.EQ.0)  GO TO 559
            MILMOT='MAILLAGE'
         ENDIF
*
*     PUIS ON REMET DES OBJETS DANS LA PILE
*
         DO 4477 IB=1,IBDDL
            IF(IB.GE.2.AND.COEX.EQ. 1.D0) CALL ECRCHA(MOTPM(1))
            IF(IB.GE.2.AND.COEX.EQ.-1.D0) CALL ECRCHA(MOTPM(2))
            CALL ECROBJ(MILMOT,MILPOI)
            CALL ECRCHA(MOTDDL(IB))
            IF(IBDDL.EQ.1) THEN
               XRCOO=XR
            ELSE
               XRCOO=XR*XNOR(IB)
            ENDIF
            CALL ECRREE(XRCOO)
 4477    CONTINUE
         GO TO 1320
      ENDIF
C
      CALL LIROBJ('POINT   ',KPOINT,0,IRETOU)
      IF(IRETOU.NE.0) GO TO 110
      CALL LIROBJ('MAILLAGE',KOBJET,1,IRETOU)
      IF(IRETOU.EQ.0)  GO TO 559
      MELEME=KOBJET
      SEGACT MELEME
      NBSOUS=LISOUS(/1)
      IF(NBSOUS.EQ.0) GO TO 120
C     OBJET COMPLEXE
      NNR=1
      DO 130 IS=1,NBSOUS
         IPT1=LISOUS(IS)
         SEGACT IPT1
         ITYESR(**)=IPT1.ITYPEL
         NBNN      =IPT1.NUM(/1)
         NBELEM    =IPT1.NUM(/2)
         MUNESR(**)=IPT1.NUM(/2)
         IF(NNR.EQ.1)LNODSR(**)=IPT1.NUM(1,1)
         DO 140  I1=1,NBELEM
            DO 1401  I2=1,NBNN
               DO 150  I3=1,NNR
                  IF(LNODSR(I3).EQ.IPT1.NUM(I2,I1)) GO TO 1401
 150           CONTINUE
               NNR=NNR+1
               LNODSR(**)=IPT1.NUM(I2,I1)
 1401       CONTINUE
 140     CONTINUE
         SEGDES IPT1
 130  CONTINUE
      SEGDES MWREL1
      SEGDES MWREL2
      SEGDES MWREL3
      SEGDES MELEME
      GO TO 160
C
C     OBJET SIMPLE
C
 120  CONTINUE
      ITYESR(**)=ITYPEL
      NBNN      =NUM(/1)
      NBELEM    =NUM(/2)
      MUNESR(**)=NUM(/2)
      NNR=0
      DO 170  I1=1,NBELEM
         DO 1701  I2=1,NBNN
            DO 180  I3=1,NNR
               IF(LNODSR(I3).EQ.NUM(I2,I1)) GO TO 1701
 180        CONTINUE
            NNR=NNR+1
            LNODSR(**)=NUM(I2,I1)
 1701    CONTINUE
 170  CONTINUE
      SEGDES MWREL1
      SEGDES MWREL2
      SEGDES MWREL3
      SEGDES MELEME
      GO TO 160
 110  CONTINUE
C
C     OBJET POINT
C
      ITYESR(**)=1
      MUNESR(**)=1
      LNODSR(**)=KPOINT
      SEGDES MWREL1
      SEGDES MWREL2
      SEGDES MWREL3
C     FINE OPERANDO RELAZIONE
 160  CONTINUE
      ICONR=IPORE1(/1)+1
C     LETTURA + O -
      IF(ICONR.EQ.2) CALL LIRMOT(MOTPM,2,IPO,0)
      IF (IPO.EQ.0)  GO TO 300
*     LIRMOT(MOTPM,2,IPO,1) goto 559
      IF(ICONR.GT.2) CALL LIRMOT(MOTPM,2,IPO,0)
      IF (IPO.EQ.0)  GO TO 300
      COEX=1.D0
      IF (IPO.EQ.2) COEX=-1.D0
C     SI CERCA DI LEGGERE UN NUOVO OPERANDO DELLA RELAZIONE
      GO TO  200
C
C     VERIFICA CONGRUENZA OPERANDI
C
 300  CONTINUE
C
      ICONR=IPORE1(/1)
      IEL11=IPORE1(1)
      IEL12=IPORE2(1)
      IEL13=IPORE3(1)
*     on autorise maintenant le point unique qui va etre applique sur
*     chaque relation
*     old      ityes=0
*     old      munes=0
      nsor=0
      nnr=0
      do 305 io=1,iconr
*     write(ioimp,*) 'io=',io
         iel11=ipore1(io)
         iel12=ipore2(io)
         iel13=ipore3(io)
         segact iel11,iel12,iel13
         nsor1=iel11.ityesr(/1)
         nsor=max(nsor,nsor1)
         NNR1=IEL13.LNODSR(/1)
         nnr=max(nnr,nnr1)
*     old         do 315 i1=1,nsor1
         if (io.eq.1) then
            segini,iel01=iel11
            segini,iel02=iel12
         else
            do 315 i1=1,nsor1
               ityes=iel01.ityesr(i1)
               iel01.ityesr(i1)=max(iel11.ityesr(i1),ityes)
*     write(ioimp,*) 'i1=',i1,' ityes=',ityes
*     $              ,' iel11.ityesr(i1)',iel11.ityesr(i1)
               munes=iel02.munesr(i1)
               iel02.munesr(i1)=max(iel12.munesr(i1),munes)
*     write(ioimp,*) 'i1=',i1,' munes=',munes
*     $              ,' iel12.munesr(i1)',iel12.munesr(i1)
 315        continue
         endif
 305  continue
      DO 310 IO=1,ICONR
*     write(ioimp,*) 'io=',io
         IEL21=IPORE1(IO)
         IEL22=IPORE2(IO)
         IEL23=IPORE3(IO)
         NSOR2=IEL21.ITYESR(/1)
*     write(ioimp,*) 'nsor,nsor2=',nsor,nsor2
         IF(NSOR.NE.NSOR2.and.nsor2.ne.1) GO TO 556
         DO 320 I2=1,NSOR2
            ityes=iel01.ityesr(i2)
            munes=iel02.munesr(i2)
*     write(ioimp,*) 'i2=',i2,' ityes=',ityes,'munes=',munes
*     write(ioimp,*) 'iel21.ityesr=',iel21.ityesr(i2)
*     $           ,' IEL22.MUNESR(I2)=',IEL22.MUNESR(I2)
            if (iel21.ityesr(i2).eq.1.and.iel22.munesr(i2).eq.1.and.
     >           nsor2.eq.1) goto 320
            IF(IEL21.ITYESR(I2).NE.ityes) GO TO 556
            IF(IEL22.MUNESR(I2).NE.munes) GO TO 556
 320     CONTINUE
         NNR2=IEL23.LNODSR(/1)
*     write(ioimp,*) 'nnr,nnr2=',nnr,nnr2
         IF(NNR2.NE.NNR.and.nnr2.ne.1) GO TO 556
 310  CONTINUE
      if (iel01.ne.0) SEGSUP,IEL01
      if (iel02.ne.0) SEGSUP,IEL02
      NBRELA=NNR
C     OPERANDI  CONGRUENTI  FINE VERIFICA
C
C     CARICAMENTO MATRICE  RIGIDEZZA
C
      NRIGEL=1
      SEGINI MRIGID
      ICHOLE=0
      IMGEO1=0
      IMGEO2=0
      ISUPEQ=0
      IFORIG=IFOUR
      COERIG(1)=1.D0
      MTYMAT='RIGIDITE'
      KRIGI=MRIGID
C
C     ON INITIALISE LE SEGMENT MELEME ASSOCIE AUX BLOCAGES
C
      SEGACT MCOORD*MOD
      NBPTSO=nbpts
*      write(ioimp,*) 'nbptso=',nbptso
      NBPTS=NBPTSO+NBRELA
      SEGADJ MCOORD
      NBSOUS=0
      NBREF=0
      NBNN=1+ICONR
      NBELEM=NBRELA
      SEGINI IPT2
      IRIGEL(1,1)=IPT2
      IPT2.ITYPEL=22
C     SI CREANO UNO PUNTO PER OGNI RELAZIONE
      DO 400 I4=1,NBRELA
         IPT2.ICOLOR(I4)=IDCOUL
         IPTS=NBPTSO+I4
         JPTS=(IPTS-1)*(IDIM+1)
         DO 410 iidim=1,idim+1
            XCOOR(JPTS+iidim)=0.D0
 410     CONTINUE
C     PUNTI ASSOCIATI AI MOLTIPLICATORI
         IPT2.NUM(1,I4)=IPTS
 400  CONTINUE
C     CARICAMENTO NODI PSEUDO-ELEMENTI
      DO 420 I8=1,ICONR
         MWREL3=IPORE3(I8)
         I10=I8+1
         SEGACT MWREL3
         LNOMAX=LNODSR(/1)
         DO 430 I9=1,NBRELA
            NN=LNODSR(min(I9,lnomax))
            NPN=(NN-1)*(IDIM+1)
            IPT2.NUM(I10,I9)=NN
            NPL1=(IPT2.NUM(1,I9)-1)*(IDIM+1)
            do 4301 iidim=1,idim+1
               XCOOR(NPL1+iidim)=XCOOR(NPL1+iidim)+XCOOR(NPN+iidim)
 4301       continue
 430     CONTINUE
         SEGDES MWREL3
 420  CONTINUE
C
C     COORDINATE DEI BARICENTRI ASSOCIATI ALLE RELAZIONI
C
*      write(ioimp,*) 'iconr,nbrela,idim=',iconr,nbrela,idim
      DO 425 I4=1,NBRELA
         NPL1=(IPT2.NUM(1,I4)-1)*(IDIM+1)
         do 4251 iidim=1,idim+1
            XCOOR(NPL1+iidim)=XCOOR(NPL1+iidim)/ICONR
 4251    continue
 425  CONTINUE

      SEGDES IPT2
      IRIGEL(2,1)=0
      IRIGEL(5,1)=NIFOUR
      IRIGEL(6,1)=NILATE
      NLIGRP=ICONR+1
      NLIGRD=NLIGRP
      SEGINI DESCR
      IRIGEL(3,1)=DESCR
      LISINC(1)='LX'
      LISDUA(1)='FLX'
      NOELEP(1)=1
      NOELED(1)=1
      DO 700 I1=1,ICONR
         I2=I1+1
         I3=ABS(INCREL(I1))
         LISINC(I2)=NOMDD(I3)
         LISDUA(I2)=NOMDU(I3)
         NOELEP(I2)=I2
         NOELED(I2)=I2
 700  CONTINUE
      SEGDES DESCR
      NELRIG=NBRELA
      rigrel=0
      SEGINI xMATRI
      IRIGEL(4,1)=xMATRI
*     LVAL=(NBNN*NBNN+NBNN)/2
      NLIGRP=NBNN
      NLIGRD=NBNN
*     SEGINI XMATRI
      DO 740 I6=1,NELRIG
*     IMATTT(I1)=XMATRI
*     740   CONTINUE
         RE(1,1,i6)= 0.D0
         I3=3
         DO 760 I1=2,NBNN
            I4=I1-1
            I2=1
            RE(I1,I2,i6)=COEFR(I4)
            RE(I2,I1,i6)=COEFR(I4)
 760     CONTINUE
 740  continue
      SEGDES XMATRI
*     SEGDES IMATRI
      call relasi(mrigid)
      SEGDES MRIGID
      CALL ECROBJ('RIGIDITE',KRIGI)
 559  CONTINUE
      ICONR=IPORE1(/1)
      DO 558 I1=1,ICONR
         MWREL1=IPORE1(I1)
         MWREL2=IPORE2(I1)
         MWREL3=IPORE3(I1)
         SEGSUP MWREL1
         SEGSUP MWREL2
         SEGSUP MWREL3
 558  CONTINUE
      SEGSUP MWGGM1
      SEGSUP MWGGM2
      SEGSUP MWGGM3
      SEGSUP MWGGM4
      SEGSUP MWGGM5
      RETURN
 556  CONTINUE
C     SEGDES IEL11,IEL12,IEL13,IEL21,IEL22,IEL23
      CALL ERREUR(324)
      GO TO 559
*
*  on arrive en 500 avec la syntaxe chp1 DUAL chp2
*
 500  continue
      CALL LIROBJ('CHPOINT',IPOCHD,1,IRECHD)
      if (ierr.ne.0) return
* ipochp est le chpoin de la relation
* ipochd est le chpoin du vecteur de controle
* decompte des dimensions
      nbsous=0
      nbref=0
      nbelem=1
      mchpoi=ipochp
      segact mchpoi
*  reserver la place pour le lx
      nligrp=1
      nbnp=0
      do ims=1,ipchp(/1)
       msoupo=ipchp(ims)
       segact msoupo
       mpoval=ipoval
       segact mpoval
       nligrp=nligrp+vpocha(/1)*vpocha(/2)
       meleme=igeoc
       segact meleme
       nbnp=nbnp+num(/2)
      enddo
*     write(6,*) ' nligrp nbelep ',nligrp,nbelep
      mchpoi=ipochd
      segact mchpoi
      nbnd=0
      xnormd = 0.d0
*  reserver la place pour le flx
      nligrd=1
      do ims=1,ipchp(/1)
       msoupo=ipchp(ims)
       segact msoupo
       mpoval=ipoval
       segact mpoval
       nligrd=nligrd+vpocha(/1)*vpocha(/2)
       do j=1,vpocha(/2)
       do i=1,vpocha(/1)
        xnormd=max(xnormd,abs(vpocha(i,j)))
       enddo
       enddo
       if (xnormd.lt.1d-5.or.xnormd.gt.1E5) then
         reaerr(1)=xnormd
         call erreur(1118)
         return
       endif
       meleme=igeoc
       segact meleme
       nbnd=nbnd+num(/2)
      enddo
*     write(6,*) ' nligrd nbeled ',nligrd,nbeled
* creation rigidite
      nligrp=nligrp+nligrd-1
      nligrd=nligrp
      segini descr
      nbnn=nbnp+1
*  maillage primal uniquement
      nbnn=nbnd+1
*  maillage dual   uniquement
      nbnn=nbnp+nbnd+1
      nbref=0
      segini meleme
      itypel=22
      nelrig=1
      rigrel=0
      segini xmatri
      symre=2
      nrigel=1
      segini mrigid
      ICHOLE=0
      IMGEO1=0
      IMGEO2=0
      ISUPEQ=0
      IFORIG=IFOUR
      COERIG(1)=1.D0
      MTYMAT='RIGIDITE'
      irigel(1,1)=meleme
      irigel(3,1)=descr
      irigel(4,1)=xmatri
      IRIGEL(5,1)=NIFOUR
      irigel(6,1)=nilate
      irigel(7,1)=2
*  remplissage maillage descripteur et valeur
*  le premier noeud sera fait a la fin. C'est le multiplicateur de lagrange
      iel=1
      ire=1
      ire=1
      nbnp=1
      nbnd=1
      mchpoi=ipochp
      do ims=1,ipchp(/1)
       msoupo=ipchp(ims)
       segact msoupo
       ipt1=igeoc
       segact ipt1
       mpoval=ipoval
       segact mpoval
       do i=1,ipt1.num(/2)
        iel=iel+1
        num(iel,1)=ipt1.num(1,i)
        nbnp=nbnp+1
       do j=1,nocomp(/2)
        ire=ire+1
        lisinc(ire)=nocomp(j)
        if (lisinc(ire).eq.'LX  ') call erreur(1125)
        noelep(ire)=iel
        re(1,ire,1)=vpocha(i,j)
        call place(nomdd,lnomdd,ipo,nocomp(j))
        if(ipo.eq.0) then
        moterr(1:4)=nocomp(j)
        moterr(5:11)='PRIMALE'
        call erreur(1117)
        endif
        lisdua(ire)=nomdu(ipo)
        noeled(ire)=iel
*       write(6,*) ' 1 ire inc dua ',lisinc(ire),lisdua(ire)
       enddo
       enddo
      enddo
      mchpoi=ipochd
      do ims=1,ipchp(/1)
       msoupo=ipchp(ims)
       segact msoupo
       ipt1=igeoc
       segact ipt1
       mpoval=ipoval
       segact mpoval
       do i=1,ipt1.num(/2)
        iel=iel+1
        num(iel,1)=ipt1.num(1,i)
        nbnd=nbnd+1
       do j=1,nocomp(/2)
        ire=ire+1
        lisdua(ire)=nocomp(j)
        noeled(ire)=iel
        re(ire,1,1)= -vpocha(i,j)
        call place(nomdu,lnomdu,ipo,nocomp(j))
        if(ipo.eq.0) then
        moterr(1:4)=nocomp(j)
        moterr(5:10)='DUALE'
        call erreur(1117)
        endif
        lisinc(ire)=nomdd(ipo)
        if (lisinc(ire).eq.'LX  ') call erreur(1125)
        noelep(ire)=iel
*       write(6,*) ' 2 ire inc dua ',lisinc(ire),lisdua(ire)
       enddo
       enddo
      enddo
*  multiplicateur de lagrange
      segact MCOORD*MOD
      nbpts=nbpts+1
      segadj mcoord
      xl=0.d0
      yl=0.d0
      zl=0.d0
      dl=0.d0
      do i=2,num(/1)
       ip=num(i,1)-1
       xl=xl+xcoor(ip*(idim+1)+1)
       if (idim.gt.1) yl=yl+xcoor(ip*(idim+1)+2)
       if (idim.gt.2) zl=zl+xcoor(ip*(idim+1)+3)
       dl=dl+xcoor(ip*(idim+1)+idim+1)
      enddo
      nbnnr=num(/1)-1
      xl=xl/nbnnr
      yl=yl/nbnnr
      zl=zl/nbnnr
      dl=dl/nbnnr
      xcoor((nbpts-1)*(idim+1)+1)=xl
      if (idim.gt.1) xcoor((nbpts-1)*(idim+1)+2)=yl
      if (idim.gt.2) xcoor((nbpts-1)*(idim+1)+3)=zl
      xcoor((nbpts-1)*(idim+1)+idim+1)=dl
      lisinc(1)= 'LX'
      lisdua(1)='FLX'
      num(1,1)=nbpts
      noelep(1)=1
      noeled(1)=1
      re(1,1,1)=0.D0
*  un resultat
      call ecrobj('RIGIDITE',mrigid)
      SEGSUP MWGGM1
      SEGSUP MWGGM2
      SEGSUP MWGGM3
      SEGSUP MWGGM4
      SEGSUP MWGGM5
      return
      END