C BLOQUE    SOURCE    PV090527  26/04/30    21:15:10     12529          

C-----------------------------------------------------------------------
C  Cet operateur impose les BLOCAGES
C
C  Syntaxe :
C
C        ENC1 = BLOQUER ( DEPL ) ( ROTA )  POI1
C
C    ou  ENC1 = BLOQUER  RADIAL  P1  (P2)  MELEME
C                        ORTHOR  P1  (P2)  MELEME
C
C    ou  ENC1 = BLOQUER (DEPL) (ROTA) DIRECTION V1  MELEME
C
C        DIM = 1        ( UX UY UZ ) ou ( UR UZ )  |
C        DIM = 2 OU 3   ( UX UY UZ RX RY RZ )      | MELEME
C        AXISYM         ( RX RZ RT UT )            |
Clist t
C        POI1   = OBJET DE TYPE POINT    (entree)
C        MELEME = OBJET DE TYPE MELEME   (entree)
C        ENC1   = OBJET DE TYPE RIGIDITE (sortie)
C
C  1) On peut imposer des BLOQUAGES UNILATERAUX en specifiant les
C     mots-cles MINIMUM ou MAXIMUM.
C  2) On peut imposer des BLOQUAGES de type FROTTEMENT en specifiant
C     le mot-cle FROTTEMENT.
C-----------------------------------------------------------------------
C  Juillet 2003 : passage a un seul multiplicateur
C-----------------------------------------------------------------------

      SUBROUTINE BLOQUE

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


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

-INC SMCHPOI
-INC TMTRAV
-INC SMCOORD
-INC SMELEME
-INC SMLMOTS
      POINTEUR MDIRCP.MLMOTS
      POINTEUR MDIRCD.MLMOTS
-INC SMMODEL
-INC SMRIGID
-INC SMTABLE

      SEGMENT MSWBLO
      CHARACTER*(LOCHPO) MOTDDL(0)
      ENDSEGMENT

      segment lispoi
        INTEGER pilpoi(mpoin),pilmul(mpoin)
      endsegment

      DIMENSION XNOR(3),U1(3),U2(3)

      CHARACTER*(LOCHPO) CHADDL,BLAN
      CHARACTER*8 MOTRIG
      CHARACTER*4 MOTPV(3), MOTBLO(5)
      CHARACTER*4 MODEPL(6),MODEDU(6), MOROTA(5),MORODU(5)
      CHARACTER*4 MODE1D(2),MOFO1D(2)

      DATA EPSI / 1.D-12 /
      DATA MOTRIG / 'RIGIDITE' /
      DATA BLAN   / '        ' /
      DATA MOTPV / 'MINI','MAXI','FROT' /
      DATA MOTBLO / 'DEPL','ROTA','RADI','ORTH','DIRE' /
      DATA MODEPL / 'UX  ','UY  ','UZ  ','UR  ','UZ  ','UT  ' /
      DATA MODEDU / 'FX  ','FY  ','FZ  ','FR  ','FZ  ','FT  ' /
      DATA MOROTA / 'RX  ','RY  ','RZ  ','RT  ','RS  ' /
      DATA MORODU / 'MX  ','MY  ','MZ  ','MT  ','MS  ' /
C Tableaux MODE1D et MOFO1D sont utilises pour certains modes 1D
      DATA MODE1D / 'UX  ','UZ  ' /
      DATA MOFO1D / 'FX  ','FZ  ' /

C  Pour ne pas avoir de verrouillage sur MCOORD en //
      SEGDES,MCOORD
      SEGACT,MCOORD*MOD

c lecture table liaisons statiques
      CALL LIRTAB('LIAISONS_STATIQUES',ipt,0,iretou)
c traitement table liaisons statiques
      IF (iretou.NE.0) THEN
        CALL BLOQU2(ipt)
        RETURN
      ENDIF

C  Est-ce une condition unilaterale ?
      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  Pas de frottement en 1D
      IF (IPO.EQ.3.AND.IDIM.EQ.1) THEN
         INTERR(1)=IDIM
         MOTERR(1:4)=MOTPV(3)
         CALL ERREUR(971)
         RETURN
      ENDIF

      idimp1=IDIM+1
C  Quelques initialisations selon le type de probleme
C  Cas IDIM = 1 :
C  ISPE1D = 1 si IDIM=1 et IFOUR=9 ou 10, car les noms de DDL primaux et
C  variables duales ne sont pas dans l'ordre "classique" (un traitement
C  specifique est alors necessaire).
      ISPE1D=0
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=1
         IAROTA=3
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
         IF (IFOUR.LE.6) THEN
            LDEPL=1
            IADEPL=0
         ELSE IF (IFOUR.GE.7.AND.IFOUR.LE.10) THEN
            LDEPL=2
            IADEPL=0
            IF (IFOUR.EQ.9.OR.IFOUR.EQ.10) ISPE1D=1
         ELSE IF (IFOUR.EQ.11) THEN
            LDEPL=3
            IADEPL=0
         ELSE IF (IFOUR.EQ.15) THEN
            LDEPL=2
            IADEPL=3
         ELSE
            LDEPL=1
            IADEPL=3
         ENDIF
         LROTA=0
         IAROTA=0
C  Autres cas :
      ELSE
         LDEPL=0
         IADEPL=0
         LROTA=0
         IAROTA=0
      ENDIF

C  Initialisation de la liste des ddls MOTDDL (segment MSWBLO)
      SEGINI,MSWBLO

C-----------------------------------------------------------------------
C  Lecture eventuelle des MOTS autres que des DDLS
C-----------------------------------------------------------------------
C     Lecture eventuelle de 'RADI','ORTH'
      IRADIA=0
      IDIREC=0
      CALL LIRMOT(MOTBLO(3),2,IMOT,0)
      IF (IMOT.NE.0) THEN
C  En DIMENSION 1, les mots-cles 'RADI,'ORTH' sont interdits
         IF (IDIM.EQ.1) THEN
            INTERR(1)=IDIM
            MOTERR(1:4)=MOTBLO(2+IMOT)
            CALL ERREUR(971)
            GOTO 1000
         ENDIF
         GOTO 44803
      ENDIF
C     Lecture eventuelle de 'DEPL' et/ou 'ROTA'
C --------------------
      IDEPL=0
      IROTA=0
 4480 CONTINUE
      CALL LIRMOT(MOTBLO,2,IMOT,0)
      IF (IMOT.EQ.1) THEN
         IDEPL=1
C  Cas particulier pour certains modes de IDIM=1
         IF (ISPE1D.EQ.1) THEN
            DO i=1,LDEPL
               MOTDDL(**)=MODE1D(IADEPL+i)
               MOTDDL(**)=MOFO1D(IADEPL+i)
            ENDDO
C  Cas general
         ELSE
            DO i=1,LDEPL
               MOTDDL(**)=MODEPL(IADEPL+i)
               MOTDDL(**)=MODEDU(IADEPL+i)
            ENDDO
         ENDIF
         GOTO 4480
      ELSEIF (IMOT.EQ.2) THEN
         IROTA=1
         DO i=1,LROTA
            MOTDDL(**)=MOROTA(IAROTA+i)
            MOTDDL(**)=MORODU(IAROTA+i)
         ENDDO
         GOTO 4480
      ENDIF
C  Lecture eventuelle de 'DIRE'
C     --------------------
      IDIREC=0
      CALL LIRMOT(MOTBLO(5),1,IMOT,0)
      IF (IMOT.NE.0) THEN
         GOTO 44804
      ENDIF
      IBDDL=MOTDDL(/2)
      IF (IBDDL.NE.0) GOTO 449
      IF (IBDDL.EQ.0) GOTO 445
C On a trouve le mot RADial ou le mot ORTHoradial
44803 CONTINUE
      IRADIA=IMOT
      CALL LIROBJ('POINT',KPOINT,1,IRETOU)
      IF (IRETOU.EQ.0) GOTO 1000
      j=(KPOINT-1)*idimp1
      DO i=1,IDIM
         U1(i)=XCOOR(j+i)
      ENDDO
      IF (IDIM.EQ.3) THEN
C  Lecture du 2nd point de l'axe (en 3D)
         CALL LIROBJ('POINT',KPOINT,1,IRETOU)
         IF (IRETOU.EQ.0) GOTO 1000
         j=(KPOINT-1)*idimp1
         YL=0.D0
         DO i=1,IDIM
            U2(i)=XCOOR(j+i)-U1(i)
            YL=YL+U2(i)*U2(i)
         ENDDO
C  Calcul du vecteur directeur unitaire de l'axe (U2)
         IF (YL.LT.EPSI) THEN
            CALL ERREUR(237)
            GOTO 1000
         ENDIF
         YL=1.D0/SQRT(YL)
         DO i=1,IDIM
            U2(i)=U2(i)*YL
         ENDDO
      ENDIF
      GOTO 449
C     On a trouve le mot DIREction
44804 CONTINUE
      IDIREC=1
      IF (IDEPL.NE.0.AND.IROTA.NE.0) THEN
         MOTERR(1:128)='DEPL ROTA'
         CALL ERREUR(-385)
         MOTERR(1:8)='DIRE'
         CALL ERREUR(803)
         RETURN
      ENDIF
*
      CALL LIROBJ('POINT',KPOINT,0,IRETOU)
      IF (IRETOU.EQ.0) THEN
*     WRITE(IOIMP,*) 'Lecture dun chpoint'
         CALL LIROBJ('CHPOINT ',MCHPOI,1,IRETOU)
         CALL ACTOBJ('CHPOINT ',MCHPOI,1)
         IF (IERR.NE.0) GOTO 1000
         IF (IDEPL.EQ.0.AND.IROTA.EQ.0) THEN
* Les composantes seront donnees par le champ
            IDICHP=2
            CALL EXTR11(MCHPOI,MDIRCD)
            SEGACT MDIRCD
            NC=MDIRCD.MOTS(/2)
            IF (NC.LE.0) THEN
               MOTERR(1:8)='CHPOINT'
* Une donnee de type %M1:8 est vide
               CALL ERREUR(1027)
               RETURN
            ENDIF
            SEGINI,MDIRCP=MDIRCD
            DO IC=1,NC
               IPRI=0
               IDUA=0
               CALL PLACE(NOMDU,LNOMDU,IDUA,MDIRCP.MOTS(IC))
               IF (IDUA.NE.0) THEN
                  MDIRCP.MOTS(IC)=NOMDD(IC)
               ELSE
                  CALL PLACE(NOMDD,LNOMDD,IPRI,MDIRCP.MOTS(IC))
               ENDIF
               MOTDDL(**)=MDIRCP.MOTS(IC)
               IF (IDUA.NE.0) THEN
                  MOTDDL(**)=NOMDU(IDUA)
               ELSE
                  IF (IPRI.NE.0) THEN
                     MOTDDL(**)=NOMDU(IPRI)
                  ELSE
                     MOTDDL(**)=BLAN
                  ENDIF
               ENDIF
            ENDDO
         ELSE
            IDICHP=1
            JGN=LOCHPO
            IBDDL=MOTDDL(/2)
            JGM=IBDDL/2
            SEGINI,MDIRCP,MDIRCD
            DO IGM=1,JGM
               IDX=((IGM-1)*2)+1
               MDIRCP.MOTS(IGM)=MOTDDL(IDX)
               MDIRCD.MOTS(IGM)=MOTDDL(IDX+1)
            ENDDO
         ENDIF
         CALL NOMC2(MCHPOI,MDIRCD,MDIRCP,MCHPO1)
         IF (IERR.NE.0) RETURN
      ELSE
         j=(KPOINT-1)*idimp1
         YL=0.D0
         DO i=1,IDIM
            XNOR(i)=XCOOR(j+i)
            YL=YL+XNOR(i)*XNOR(i)
         ENDDO
         IF (YL.LT.EPSI) THEN
            CALL ERREUR(239)
            GOTO 1000
         ENDIF
         YL=1.D0/SQRT(YL)
         DO i=1,IDIM
            XNOR(i)=XNOR(i)*YL
         ENDDO
         IDICHP=0
      ENDIF
      GOTO 449
C
C  Lecture eventuelle de DDLs :
C ------------------------------
C
 445  CONTINUE
C  La liste des grandeurs PRIMALEs/DDLs est un LISTMOTS.
cbp,2020-11-12 : possibilite non documentee => ajout a la notice
c                + ajout de qq lignes pour rendre listmot dual facultatif
      CALL LIROBJ('LISTMOTS',MLMOT1,0,IRETOU)
      IF (IERR.NE.0) GOTO 1000
      IF (IRETOU.EQ.0) GOTO 446
      CALL LIROBJ('LISTMOTS',MLMOT2,0,IRETO2)
      IF (IERR.NE.0) GOTO 1000
C    -cas ou la liste des grandeurs DUALEs n'est pas fournie.
      IF(IRETO2.EQ.0) THEN
        SEGACT,MLMOT1
        IMOT1 = MLMOT1.MOTS(/2)
        IF (IMOT1.LE.0) THEN
          CALL ERREUR(643)
          SEGDES,MLMOT1
          GOTO 1000
        ENDIF
        DO IMOT = 1, IMOT1
          MOTDDL(**)=MLMOT1.MOTS(IMOT)
          CALL PLACE(NOMDD,LNOMDD,idd,MLMOT1.MOTS(IMOT))
          IF (IERR.NE.0) GOTO 1000
          IF (idd.LE.0) THEN
             MOTDDL(**)=BLAN
          ELSE
             MOTDDL(**)=NOMDU(idd)
          ENDIF
        ENDDO
        SEGDES,MLMOT1
C    -cas ou la liste des grandeurs DUALEs est un LISTMOTS fourni.
cbp,2020-11-12 : on ne documente pas cette possibilite car pas le temps
      ELSE
        SEGACT,MLMOT1,MLMOT2
        IMOT1 = MLMOT1.MOTS(/2)
        IF (IMOT1.LE.0) THEN
          CALL ERREUR(643)
          SEGDES,MLMOT1,MLMOT2
          GOTO 1000
        ELSE IF (IMOT1.NE.MLMOT2.MOTS(/2)) THEN
          CALL ERREUR(854)
          SEGDES,MLMOT1,MLMOT2
          GOTO 1000
        ENDIF
        DO IMOT = 1, IMOT1
          MOTDDL(**)=MLMOT1.MOTS(IMOT)
          MOTDDL(**)=MLMOT2.MOTS(IMOT)
        ENDDO
        SEGDES,MLMOT1,MLMOT2
      ENDIF
      GOTO 449

 446  CONTINUE
C  La liste des grandeurs PRIMALEs/DDLs sont des MOTs.
C  La liste des ddls autorises = NOMDD defini dans BDATA.ESO
C  On doit lire au moins 1 ddl (car sinon MOTDDL est vide !)
      CALL LIRMOT(NOMDD,LNOMDD,IMOT,0)
      IF (IERR.NE.0) GOTO 1000
      IF (IMOT.NE.0) THEN
         MOTDDL(**)=NOMDD(IMOT)
         MOTDDL(**)=NOMDU(IMOT)
      ELSE
         CALL LIRCHA(CHADDL,0,IMOT)
         IF (IERR.NE.0) GOTO 1000
         IF (IMOT.EQ.0) GOTO 447
         MOTDDL(**)=CHADDL
         MOTDDL(**)=BLAN
      ENDIF
      GOTO 446

C  Conditions aux limites via un MODELE --> non documente, non plus !
 447  CONTINUE
      IBDDL=MOTDDL(/2)
      IPMODL = 0
      CALL LIROBJ('MMODEL  ',IPMODL,0,IRETOU)
      IF (IRETOU.EQ.0) GOTO 449
      CALL ACTOBJ('MMODEL  ',IPMODL,1)
      CALL NOVARD(IPMODL,'DEPL')
      CALL LIROBJ('LISTMOTS',MLMOT1,1,IRETOU)
      IF (IERR.NE.0) GOTO 1000
      CALL NOVARD(IPMODL,'FORC')
      CALL LIROBJ('LISTMOTS',MLMOT2,1,IRETOU)
      IF (IERR.NE.0) GOTO 1000
      SEGACT,MLMOT1,MLMOT2
      IMOT1 = MLMOT1.MOTS(/2)
      DO i = 2, IBDDL, 2
         IF (MOTDDL(i).EQ.BLAN) THEN
            CALL PLACE(MLMOT1.MOTS(1),IMOT1,iplac,MOTDDL(i-1))
            IF (iplac.EQ.0) THEN
               MOTERR(1:4)=MOTDDL(i-1)
               CALL ERREUR(197)
               IF (IERR.NE.0) GOTO 1000
            ELSE
               MOTDDL(i)=MLMOT2.MOTS(iplac)
            ENDIF
         ENDIF
      ENDDO
      SEGSUP,MLMOT1,MLMOT2
 449  CONTINUE
      IBDDL=MOTDDL(/2)
      DO i = 2, IBDDL, 2
         IF (MOTDDL(i).EQ.BLAN) THEN
            MOTERR(1:4)=MOTDDL(i-1)
            CALL ERREUR(108)
            IF (IERR.NE.0) GOTO 1000
         ENDIF
      ENDDO
C  Verification que le nombre de DDLs a bloquer n'est pas nul
C     IF (IBDDL.EQ.0) GOTO 1000
C-----------------------------------------------------------------------
C  Fin de la lecture des mots (DEPL,ROTA...) ou des DDLs
C-----------------------------------------------------------------------

C  Recherche du maillage MELEME de type POINT :
C ----------------------------------------------
      KOBJET=0
      MELEME=0
C  On cherche d'abord si on a un POINT que l'on transformera en POI1
C  sinon on cherche un maillage que l'on transforme en POI1 si besoin
      CALL LIROBJ('POINT',KPOINT,0,IRETOU)
      IF (IERR.NE.0) GOTO 1000
      IF (IRETOU.EQ.0) THEN
         CALL LIROBJ('MAILLAGE',KOBJET,1,IRETOU)
         IF (IERR.NE.0) GOTO 1000
         MELEME=KOBJET
         SEGACT,MELEME
         IF (ITYPEL.NE.1) CALL CHANGE(MELEME,1)
         NBPOIN=NUM(/2)
      ELSE
         MELEME=KPOINT
         CALL CRELEM(MELEME)
         NBPOIN=1
      ENDIF
*
      IF (IDICHP.GE.1) THEN
* Cas ou la direction est donnee par un chpoint
* On construit un segment MTRAV contenant les directions normees
* (DEPL ou ROTA) ou non
*         WRITE(IOIMP,*) 'Transfo du chpoint en MTRA'
* 1) Recopie des composantes et valeurs pertinentes du chpoint
*     dans le TMTRAV
         CALL CP2TR2(MDIRCP,MELEME,MCHPO1,MTRAV)
         IF (IERR.NE.0) GOTO 1000
         SEGACT MELEME
         SEGACT MTRAV
*         segprt,MTRAV
* 2) On norme la direction dans le cas DEPL ou ROTA
*         WRITE(IOIMP,*) 'Calcul des normes '
         IF (IDICHP.EQ.1) THEN
            SEGACT MTRAV*MOD
            DO IBPOIN=1,NBPOIN
               YL=0.D0
               DO I=1,IDIM
                  XNOR(I)=BB(I,IBPOIN)
                  YL=YL+XNOR(I)*XNOR(I)
               ENDDO
               IF (YL.LT.XPETIT) THEN
                  CALL ERREUR(239)
                  GOTO 1000
               ENDIF
               YL=1.D0/SQRT(YL)
               DO I=1,IDIM
                  BB(I,IBPOIN)=XNOR(I)*YL
               ENDDO
            ENDDO
         ENDIF
*         WRITE(IOIMP,*) 'Fin du calcul des normes '
*         SEGPRT,MTRAV
      ENDIF

C  Determination du nombre de multiplicateurs NNMAT par noeud de MELEME
C  NNMAT correspond au nombre de DDLs a bloquer par noeud = nombre de
C  multiplicateurs a creer par noeud (1 multiplicateur) = NRIGEL
C  Dans les cas RADIal et DIREction, on a une seule matrice par noeud.
C  Dans les autres cas, autant de matrices que MOTDDL(/2)/2.
      IF (IDIREC.NE.0.OR.IRADIA.NE.0) THEN
         NNMAT=1
      ELSE
         NNMAT=IBDDL/2
      ENDIF
C  Initialisation de l'objet RIGIDITE associe aux BLOCAGES
      NRIGEL=NNMAT
      SEGINI,MRIGID
      MTYMAT=MOTRIG
      IFORIG=IFOUR
      ICHOLE=0
      IMGEO1=0
      IMGEO2=0
C  A chaque multplicateur est associe un nouveau noeud
C  NBPOIN : nombre de points du maillage MELEME a bloquer
      NBNO=nbpts
      NBNOI=NBNO
      NBPTS=NBNO+NNMAT*NBPOIN
      SEGADJ,MCOORD

C Boucle sur le nombre de DDLs a bloquer
      DO IAA=1,NNMAT
C  Pour chaque noeud du maillage, creation d'un noeud associe (place
C  initialement a l'origine) au IAA-eme DDL
         DO i=1,NBPOIN
            DO j=1,idimp1
               XCOOR(NBNOI*idimp1+j)=0.D0
            ENDDO
            NBNOI=NBNOI+1
         ENDDO
C  Creation du maillage MELEME de MULTiplicateurs associe aux BLOCAGES
         NBSOUS=0
         NBREF=0
         NBNN=2
         NBELEM=NBPOIN
         SEGINI,IPT1
         KIPT1=IPT1
         IPT1.ITYPEL=22
         DO i=1,NBPOIN
            j=(IAA-1)*NBPOIN+i
            IPT1.NUM(1,i)=NBNO+j
            IPT1.NUM(2,i)=NUM(1,i)
            IPT1.ICOLOR(i)=IDCOUL
C  Correction : Les pts mult sont a la meme position que les noeuds
            IREF3=(NUM(1,i)-1)*idimp1
            IREF1=(NBNO+j-1)*idimp1
            DO j=1,IDIM
               XCOOR(IREF1+j)=XCOOR(IREF3+j)
            ENDDO
         ENDDO
C  Creation des RAIDEURS associees au IAA-eme multplicateur (DDL)
         IRIGEL(1,IAA)=KIPT1
         IRIGEL(2,IAA)=0
         IRIGEL(5,IAA)=NIFOUR
         IRIGEL(6,IAA)=NILATE
         IF (IDIREC.NE.0.OR.IRADIA.NE.0) THEN
            IF (IDICHP.EQ.2) THEN
               NCBLO=MDIRCP.MOTS(/2)
            ELSE
               NCBLO=LDEPL
            ENDIF
         ELSE
            NCBLO=1
         ENDIF
         NLIGRP=1+NCBLO
         NLIGRD=NLIGRP
C  Remplissage du tableau des DESCripteurs de RIG
         SEGINI,DESCR
         IRIGEL(3,IAA)=DESCR
         NOELEP(1)=1
         NOELEP(2)=2
         NOELED(1)=1
         NOELED(2)=2
         IF (IDIREC.NE.0.OR.IRADIA.NE.0) THEN
            IF (IDICHP.EQ.2) THEN
               DO i=1,NCBLO
                  NOELEP(1+i)=2
                  NOELED(1+i)=2
                  LISINC(1+i)=MOTDDL(2*i-1)
                  LISDUA(1+i)=MOTDDL(2*i)
               ENDDO
            ELSE
               DO i=1,NCBLO
                  NOELEP(1+i)=2
                  NOELED(1+i)=2
                  IF (IROTA.NE.1) THEN
                     LISINC(1+i)=MODEPL(IADEPL+i)
                     LISDUA(1+i)=MODEDU(IADEPL+i)
                  ELSE
                     LISINC(1+i)=MOROTA(IADEPL+i)
                     LISDUA(1+i)=MORODU(IADEPL+i)
                  ENDIF
               ENDDO
            ENDIF
         ELSE
            j=2*(IAA-1)
            LISINC(2)=MOTDDL(j+1)
            LISDUA(2)=MOTDDL(j+2)
         ENDIF
         LISINC(1)='LX  '
         LISDUA(1)='FLX '
*         SEGDES DESCR
         NELRIG=NBPOIN
         rigrel=0
         SEGINI,XMATRI
         IRIGEL(4,IAA)=XMATRI
         COERIG(IAA)=1.D0
C  Remplissage de  la matrice de rigidite RE :
C  Il faut distinguer les cas IRADIA, IDIREC et les autres
C  IRADIA : Il faut calculer la DIREction puis identique a IDIREC
C  IDIREC : La DIRECTION est stockee dans le vecteur norme XNOR
C  AUTRES : La matrice est predefinie dans RIG
C  Option RADIAL : Calcul prealable de la direction pour chaque noeud
         IF (IRADIA.NE.0) THEN
            DO IB=1,NBPOIN
               j=(NUM(1,IB)-1)*idimp1
               DO i=1,IDIM
                  XNOR(i)=XCOOR(j+i)-U1(i)
               ENDDO
               IF (IDIM.EQ.2) THEN
                  YL=XNOR(1)*XNOR(1)+XNOR(2)*XNOR(2)
                  IF (YL.LT.EPSI) THEN
                     CALL ERREUR(238)
                     RETURN
                  ENDIF
                  YL=1.D0/SQRT(YL)
                  IF (IRADIA.EQ.1) THEN
                     XNOR(1)=XNOR(1)*YL
                     XNOR(2)=XNOR(2)*YL
                  ELSE IF (IRADIA.EQ.2) THEN
                     XX=XNOR(1)
                     XNOR(1)=-XNOR(2)*YL
                     XNOR(2)=XX*YL
                  ENDIF
               ELSE
                  YL=XNOR(1)*U2(1)+XNOR(2)*U2(2)+XNOR(3)*U2(3)
                  XL=0.D0
                  DO i=1,3
                     XNOR(i)=XNOR(i)-YL*U2(i)
                     XL=XL+XNOR(i)*XNOR(i)
                  ENDDO
                  IF (XL.LT.EPSI) THEN
                     CALL ERREUR(238)
                     RETURN
                  ENDIF
                  IF (IRADIA.EQ.1) THEN
                     XL=1.D0/SQRT(XL)
                     XNOR(1)=XNOR(1)*XL
                     XNOR(2)=XNOR(2)*XL
                     XNOR(3)=XNOR(3)*XL
                  ELSE IF (IRADIA.EQ.2) THEN
                     XX=XNOR(1)
                     YY=XNOR(2)
                     ZZ=XNOR(3)
                     XNOR(1)=YY*U2(3)-ZZ*U2(2)
                     XNOR(2)=ZZ*U2(1)-XX*U2(3)
                     XNOR(3)=XX*U2(2)-YY*U2(1)
                  ENDIF
               ENDIF
C  XNOR contient la direction normee
               RE(1,1,IB)=0.D0
               DO IC=1,IDIM
                  RE(IC+1,1,IB)=XNOR(IC)
                  RE(1,IC+1,IB)=XNOR(IC)
               ENDDO
            ENDDO
C  Option DIRECTION
         ELSE IF (IDIREC.EQ.1) THEN
            IF (IDICHP.GE.1) THEN
               DO I=1,NBPOIN
                  RE(1,1,i)=0.D0
                  DO IC=1,NCBLO
                     RE(IC+1,1,i)=BB(IC,I)
                     RE(1,IC+1,i)=BB(IC,I)
                  ENDDO
               ENDDO
            ELSE
               DO I=1,NBPOIN
                  RE(1,1,i)=0.D0
                  DO IC=1,IDIM
                     RE(IC+1,1,i)=XNOR(IC)
                     RE(1,IC+1,i)=XNOR(IC)
                  ENDDO
               ENDDO
            ENDIF
C  Autres options :
         ELSE
            DO i=1,NBPOIN
               RE(1,1,i)=0.D0
               RE(2,1,i)=1.D0
               RE(2,2,i)=0.D0
               RE(1,2,i)=RE(2,1,i)
            ENDDO
         ENDIF
      ENDDO
C  Fin de la boucle sur les IAA DDLs a bloquer
      IF (IDICHP.GE.1) THEN
*     WRITE(IOIMP,*) 'Destruction du segment de travail'
         SEGSUP,MTRAV
         SEGSUP MCHPO1
         SEGSUP MDIRCP
         SEGSUP MDIRCD
      ENDIF
C  Fin normale du traitement
      call relasi(mrigid)
      CALL ACTOBJ('RIGIDITE',MRIGID,1)
      CALL ECROBJ('RIGIDITE',MRIGID)
      IF (MELEME.GT.0) THEN
         IF (KOBJET.EQ.0) THEN
            SEGSUP,MELEME
         ELSE
*            SEGDES,MELEME
         ENDIF
      ENDIF
C  Sortie du sousprogramme (normale ou en cas d'erreur)
 1000 CONTINUE
      SEGSUP,MSWBLO
*      SEGDES,MCOORD
      END
 
 
 
