Sources Esope | Cast3M

assem5
C ASSEM5    SOURCE    PV        22/10/12    21:15:01     11474          
      SUBROUTINE ASSEM5(ITRAV1,ITOPO1,INUIN1,MMMTRI,INCTR1
     &,IITOP1,NBNNMA,SNTT)
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c  ****  subroutine pour faire l'assemblage de matrices symetriques
c        en vue d'une decomposition de crout modifiee.
c
c  en entree:
c  ****  itrav1 : pointeur objet mrigidite
c  ****  itopo1 : pointeur segment de travail itopo ( voir assem1)
c  ****  iitop1 : pointeur segment de travail iitop ( voir assem1)
c  ****  inuin1 : pointeur segment de travail inuinv( voir assem1)
c  ****  mmmtri : pointeur objet matrice triangularisee (non modifie)
c                 (voir smmatri)
c  ****  nbnnma : nombre d'inconnues maitresses
c  ****  sntt   : segment des noeuds non totalement maitres
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
-INC PPARAM
-INC CCOPTIO
-INC SMELEME
-INC SMRIGID
-INC SMMATRI
-INC CCREEL
      SEGMENT,INUINV(NNGLOB)
      SEGMENT,ITOPO(IENNO)
      SEGMENT,IITOP(NNOE+1)
      SEGMENT,INCTRR(NIRI)
      SEGMENT,INCTRA(NLIGRE)
      SEGMENT,IPV(NNOE)
      SEGMENT,VMAX(INC)
      SEGMENT SNTT
      INTEGER NTTMAI(NN)
      ENDSEGMENT
c
c  ****  ces tableaux servent au reperage de la matrice pour l'assemblag
c  ****  il seront tous supprimes en fin d'assemblage.
c
      SEGMENT,IVAL(NNN)
      SEGMENT,ITRA(NNN,2)
      SEGMENT TRATRA
      REAL*8  XTRA(INCRED,INCDIF)
      INTEGER LTRA(INC,INCDIF)
      INTEGER NTRA(INCRED,INCDIF)
      INTEGER MTRA(INCDIF)
      ENDSEGMENT
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c  ****  ival(i)=j     : la i eme ligne d'une petite matrice s'assemble
c                      dans la j eme de la grande.
c  ****   itrav(i,1)=j : la ieme inconnue du noeud en cours d'assemblage
c                      et qui se trouve dans la petite matrice se trouve
c                      en j eme position de la petite matrice.
c  ****   itrav(i,2)   : la ieme inconnue du noeud en cours d'assemblage
c                      present dans la petite matrice est en jeme
c                      position dans la grande
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      SEGMENT,RA(N1,N1)*D
      SEGMENT JNOMUL
      LOGICAL INOMUL(NNR)
      ENDSEGMENT
      REAL*8 DMAX,COER,DDDD,DMAXY
      LOGICAL NOMUL
*
* en cas de normalisation des variables.
*
      INWUIT=0
      IF(NORINC.NE.0) CALL ASSEM0(ITRAV1,1,INWUIT)
*
*  pv  on active une fois pour toutes les meleme descr... de la rigidite
*      on en profite pour creer inomul
c
c  ****  recherche de la dimension max de ival,et segini de ival et itra
c
c
      INCTRR=INCTR1
      SEGACT,INCTRR
c
      MRIGID=ITRAV1
 
      SEGACT,MRIGID
      NNR=IRIGEL(/2)
      NNN=0
      SEGINI JNOMUL
      DO 1 IRI=1,NNR
c
          INCTRA=INCTRR(IRI)
          SEGACT INCTRA
c
          DESCR=IRIGEL(3,IRI)
          SEGACT, DESCR
c
          IPT1=IRIGEL(1,IRI)
          SEGACT IPT1
c
          XMATRI=IRIGEL(4,IRI)
          SEGACT XMATRI
c
c    ... na = nb de variables primales ...
          NA=LISINC(/2)
c    ... nnn = max des nb de variables primales ...
          NNN=MAX(NA,NNN)
c
c    ... inomul (non multiplicateurs ?) dit s'il ne s'agit pas des cl ...
          INOMUL(IRI)=.TRUE.
          IF(IPT1.ITYPEL.EQ.49) INOMUL(IRI)=.FALSE.
c
   1  CONTINUE
c
c ... ival a la taille nnn ...
      SEGINI,IVAL
c ... itra a la taille (nnn,2) ...
      SEGINI,ITRA
c
c  ****  activation des segments de travails et de mmatri
c
c ... recherche du nombre d'inconnues impliquees
c
      MMATRI=MMMTRI
      SEGACT,MMATRI*MOD
c
      MINCPO=IINCPO
      SEGACT,MINCPO
c
c
      MAXIPO=0
      DO  I=1,INCPO(/2)
          DO  J=1,INCPO(/1)
              MAXIPO=MAX(INCPO(J,I),MAXIPO)
          ENDDO
      ENDDO
c
      N1=MAXIPO
c
      ITOPO=ITOPO1
      SEGACT,ITOPO
c
      IITOP=IITOP1
      SEGACT,IITOP
c
      INUINV=INUIN1
      SEGACT,INUINV
c
c ... nnoe = nombre de noeuds impliqués ...
c
      SEGACT SNTT
      NNOE=INCPO(/2)+NTTMAI(/1)
c
c ... inc = dimension de la matrice ...
      INC=MAXIPO
c
      SEGINI,MDIAG
      IDIAG=MDIAG
c
      SEGINI,MILIGN
      IILIGN=MILIGN
c
c ... incdif = nb de variables primales différentes ...
      INCDIF=INCPO(/1)
c
      MIMIK=IIMIK
      SEGACT,MIMIK
c
      NNOE=INCPO(/2)
c
c ... ipv est de taille nnoe ...
      SEGINI IPV
c
      INCRED=0
      DO 80 INO =1,NNOE
          ICOMPT=0
c    ... maxele = nb d'éléments qui contiennent le noeud no ino ...
          MAXELE = (IITOP(INO+1)-IITOP(INO))/2
c    ... boucle sur ces éléments ...
          DO 81 IELE=1,MAXELE
c
c       ... iiu = pointeur dans iitop ...
              IIU=IITOP(INO) + IELE + IELE -2
c
c       ... iel = numero de l'élément dans son maillage ...
              IEL=ITOPO(IIU)
c
c       ... iri = numero du maillage (dans irigel) de cet élément ...
              IRI=ITOPO(IIU+1)
c
c       ... on va chercher ce maillage dans irigel ...
              MELEME=IRIGEL(1,IRI)
c
c       ... boucle sur les noeuds de cet élément ...
              DO 83 I=1,NUM(/1)
c
c          ... ip = numéro local du i-ème noeud de l'élément ...
                  IP=INUINV(NUM(I,IEL))
c
c          ... la première condition = les "connections" sont symétriques ...
**pv              IF (IP.GT.INO) GOTO 83
c          ... le noeud ip a déjà vu ce noeud ...
                  IF (IPV(IP).EQ.INO) GOTO 83
c
c          ... ces deux opérations se font si (ip &lt;= ino et ipv(ip) != ino) ...
c          ... ipv = n° max des noeuds "connectés" au noeud ip ...
c
                  IPV(IP)=INO
c          ... icompt = nombre de noeuds connectés (de numéro &lt; ou =) au noeud ino ...
                  ICOMPT=ICOMPT+1
  83          CONTINUE
  81      CONTINUE
c    ... incred = maximum des icompt sur tous les noeuds impliqués ...
c
          INCRED=MAX(INCRED,ICOMPT)
  80  CONTINUE
c
      SEGSUP IPV
c
      INCRED=INCRED*INCDIF
c ... la taille de tratra dépend de inc, incred et incdif ...
      SEGINI TRATRA
      segini vmax
c
c ... coefficients de normalisation ...
c     -----------------------------
c
      LLVNUL=0
      IJMAX=0
      NJTOT=0
c
c ... la taille de mdnor est inc (nombre d'équations) ...
      SEGINI MDNOR
      IDNORM=MDNOR
      DO 20 I=1,INC
          DNOR(I)=1.D0
   20 CONTINUE
c
c ... si une normalisation est imposée, les valeurs dans mdnor seront modifiées
c    par assem3 ...
      IF(NORINC.NE.0) CALL ASSEM3(MDNOR,MIMIK,MINCPO)
c
c  ****  boucle   *100* sur les numeros de noeuds que l'on assemble
c        -------------------------------------------------------------
c
c    ... ipre = n° de la première ligne concernée par le noeud ino ...
      IPRE=1
c
      DO 999 IJK=1,2
c
          DO 100 INO=1,NNOE
              DO 101 IIT=1,INCDIF
                  MTRA(IIT)=0
 101          CONTINUE
c    ... ider = n° de la dernière ligne concernée par le noeud ino ...
              IDER=0
c
              IF(IJK.EQ.2) IPRE=MAXIPO+1
              DO 998 II=1,INCDIF
                  IF (IJK.EQ.1.AND.INCPO(II,INO).LE.NBNNMA)
     #            IDER=MAX(INCPO(II,INO),IDER)
                  IF (IJK.EQ.2.AND.INCPO(II,INO).GT.NBNNMA) THEN
                      IDER=MAX(INCPO(II,INO),IDER)
                      IPRE=MIN(INCPO(II,INO),IPRE)
                  ENDIF
 998          CONTINUE
              IF (IDER.EQ.0.OR.IPRE.EQ.(MAXIPO+1)) GOTO 100
              LLVVA=0
c
c  ****  boucle   *99*   sur les elements touchant le noeud ino
c                        pour les elements multiplicateur on ne fait pas
c                        l'assemblage
c
c    ... maxele = nb d'éléments qui contiennent le noeud no ino ...
              MAXELE= (IITOP(INO+1) -IITOP(INO))/2
c
 
              DO 99 IELE=1,MAXELE
                  IIU=IITOP(INO) + IELE + IELE - 2
c       ... iel = numero de l'élément dans son maillage ...
                  IEL=ITOPO(IIU)
c       ... iri = numero du maillage (dans irigel) de cet élément ...
                  IRI=ITOPO(IIU+1)
c       ... meleme = pointeur vers ce maillage ...
                  MELEME=IRIGEL(1,IRI)
c       ... descr = pointeur vers le segment descr concerné ...
                  DESCR=IRIGEL(3,IRI)
c       ... inctra contient les indices (dans imik et ihar) des descriptions
c           des ddl correspondants ...
                  INCTRA=INCTRR(IRI)
c       ... xmatri = tableau de pointeurs vers les matrices élémentaires ...
                  XMATRI=IRIGEL(4,IRI)
 
c       ... coer = coefficient multiplicateur ...
                  COER=COERIG(IRI)
c
c  ****  nomul =.false.  il existe  un multiplicateur
c  ****  initialisation de ival. ival(i)=j veut dire que
c  ****  la i eme ligne de la petite matrice s'assemble dans
c  ****  la j eme de la grande matrice.
c
c       ... nin = nombre de variables primales (dans descr) ...
                  NIN=LISINC(/2)
c       ... non multiplicateur ? ...
                  NOMUL=INOMUL(IRI)
                  NA=0
c       ... boucle sur les variables primales de l'élément ...
c
                  DO 98 ICO=1,NIN
c          ... ija = numéro local du noeud-support du ddl ...
                      IJA=INUINV(NUM(NOELEP(ICO),IEL))
c          ... ijb = indice du ddl ...
                      IJB=INCTRA(ICO)
c          ... on met dans ival(ico) le n° de l'équation correspondant au noeud ija
c              et ddl no ijb ...
                      IVAL(ICO)=INCPO(IJB,IJA)
c          ... si on a trouvé le bon noeud ...
                      IF(IJA.EQ.INO) THEN
c             ... on incrémente na (le nombre de ddl connectés au noeud ino ?) ...
                          NA=NA+1
                          ITRA(NA,1)=ICO
                          ITRA(NA,2)=IVAL(ICO)
                      ENDIF
  98              CONTINUE
*                  XMATRI=IMATTT(IEL)
*                  SEGACT,XMATRI
c
c  ****  boucle  *95* sur les inconnues de la petite matrice
c
                              DO 90 INCC=1,NA
c             ... inco = le n° d'équation (ou de ligne) du ddl no incc du noeud ino ...
                          INCO=ITRA(INCC,2)
                          IF (IJK.EQ.1) then
                              IF (INCO.GT.IDER) GOTO 90
                          ENDIF
c             ... iloc = n° de ligne relatif / ipre du noeud ino ...
c
                          ILOC=INCO-IPRE+1
                          IF (IJK.EQ.2) THEN
                              IF (ILOC.LE.0) GOTO 90
                          ENDIF
c
c             ... jj = n° (dans l'élément) de la colonne (et de la ligne) correspondant au
c                 ddl n° incc du noeud ino ...
                          JJ=ITRA(INCC,1)
                  DO 95 IK=1,NIN
c          ... io = n° de l'équation correspondant au ddl no ik de l'élément (n° de la colonne) ...
                      IO=IVAL(IK)
c          ... s'il dépasse le dernier intéressant on ne fait rien ...
                      IF(IO.GT.IDER) GO TO 95
c          ... boucle sur les ddl's du noeud ino ...
c
c
c             ... si io > inco on ne fait plus rien dans cette boucle ...
c             ... car c'est le cas symétrique !!! ...
                          IF (IO.GT.INCO)  GO TO 95
c
c             ... ci-dessous on stockera les quelques lignes de la rigidité concernant le
c                 noeud ino dans le segment tratra. les coefficients vont dans xtra, le deuxième
c                 indice (iloc) est le numéro relatif de la ligne, le premier (imtt) n'a rien à voir avec le
c                 numéro de la colonne. celui-ci est stocké dans ntra(imtt,iloc). l'information inverse
c                 (imtt en fonction de io) se trouve dans ltra(io,iloc). on met dans mtra le nombre
c                 de termes différents assemblés par ligne ...
                          ILTT= LTRA(IO,ILOC)
                          IF(ILTT.EQ.0) THEN
c                ... llvva est initialisé à 0 au début de la boucle 100 ...
                              LLVVA=LLVVA+1
c                ... les mtra sont mis à 0 au début de la boucle 100 ...
                              IMMTT=MTRA(ILOC)+1
                              MTRA(ILOC)=IMMTT
                              XTRA(IMMTT,ILOC)=0.D0
                              NTRA(IMMTT,ILOC)=IO
                              LTRA(IO,ILOC)=IMMTT
                              ILTT=IMMTT
                          ENDIF
c
c             ... si ce n'est pas une condition aux limites, on assemble ...
                          IF(NOMUL) THEN
cfaux !!!            xtra(iltt,iloc)=xtra(iltt,iloc)+re(ik,jj)*coer
c                ... cette ligne a été remplacée par ce qui suit car pour les termes
c                    de couplage entre les ddl du même noeud, les termes au-dessus de
c                    la diagonale étaient pris. Ça marchait car les matrices étaient
c                    symétriques, comme maintenant on travaillera aussi sur des asymétriques ...
***                           IF(JJ.LT.IK) THEN
***                               IIILIG=IK
***                               IIICOL=JJ
***                           ELSE
***                               IIILIG=JJ
***                               IIICOL=IK
***                           ENDIF
c                ... attention ! dans xtra 1er indice est lié à la colonne, deuxième est le numéro
c                    de la ligne, tandis que dans re c'est l'inverse ...
 
*** on exploite la symetrie de re - XTRA comprend COER
                   XTRA(ILTT,ILOC)=XTRA(ILTT,ILOC)+RE(ik,jj,IEL)
     &             *COER
                          ENDIF
   95             CONTINUE
   90                 CONTINUE
*                  SEGDES,XMATRI
   99         CONTINUE
c
c  *** compactage des lignes, en meme temps calcul de ijmax qui sera
c  *** la dimension max d'un segment lign.
c  *** le segment associe a une ligne (segment llign) est de la forme :
c  *** immmm(na) permet de savoir si un mouvement d'ensemble sur la
c  *** ligne existe. ippo(na+1) donne la position dans xxva la 1ere
c  *** valeur de la ligne .xxva  valeur de la matrice.
c  *** linc(i) donne le numero de la colonne du ieme elem de xxva
c
c    ... na = nb de lignes concernant le noeud ino ...
              NA = IDER-IPRE+1
c    ... llvnul = somme cumulée des llvva ...
              LLVNUL=LLVNUL+LLVVA
c
              SEGINI,LLIGN
c
              NBINO=INO
              IF (IJK.EQ.2) THEN
                  DO 116 JJ=1,NTTMAI(/1)
                      IF (NTTMAI(JJ).EQ.INO) GOTO 117
 116              CONTINUE
                  GOTO 118
 117              CONTINUE
                  NBINO=NNOE+JJ
              ENDIF
 118          CONTINUE
              ILIGN(NBINO)=LLIGN
              NBA=0
c
c    ... boucle sur les ddl du noeud ino ...
              DO 120 JPA=1,NA
c       ... iiin = n° de ligne du ddl n° jpa ...
                  IIIN=IPRE+JPA -1
c
c       ... que l'on met au bon endroit dans immmm (faisant partie de llign) ...
c
                  IMMMM(JPA)=IIIN
c
c
c       ... ippo(i)+1 = adresse du début (dans xxva et linc) des données relatives au ddl n° i ...
                  IPPO(JPA)=NBA
c
c       ... boucle sur les termes dans la ligne ...
                  DO 121 IPAK = 1,MTRA(JPA)
c          ... iunpak = n° de la colonne du terme n° ipak ...
                      IUNPAK=NTRA(IPAK,JPA)
c          ... pour les termes mis dans llign on efface l'information sur la position dans xtra ...
                      LTRA(IUNPAK,JPA)=0
c          ... compteur ++ ...
                      NBA=NBA+1
c          ... le n° de la colonne va dans linc ...
                      LINC(NBA)=IUNPAK
c          ... transfert du xtra (segment tratra) vers xxva (segment llign) ...
                      XXVA(NBA)=XTRA(IPAK,JPA)
             vmax(iiin)=max(abs(xxva(nba)),vmax(iiin))
c          ... les termes diagonaux vont dans diag (segment mdiag) ...
                      IF(IIIN.EQ.IUNPAK) DIAG(IIIN)=XXVA(NBA)
  121             CONTINUE
  120         CONTINUE
c    ... le pointeur vers la fin de la dernière ligne ...
              IPPO(NA+1)= NBA
              NJMAX= 0
*  recherche du mini globale sur toutes les inconnues
      LPA=IPRE
c    ... boucle sur les lignes relatives au noeud ino ...
              DO 126 JPA=IPRE,IDER
c       ... on met le n° du noeud dans ipno (segment milign) ...
                  IPNO(JPA)=NBINO
c       ... ipde et ipdf : début et fin des données relatives au ddl n° jpa dans xxva et linc ...
                  IPDE=IPPO(JPA-IPRE+1)+1
                  IPDF=IPPO(JPA-IPRE+2)
c       ... lpa = n° mini de la colonne avec des termes non nuls dans la ligne n° jpa ...
                  DO 155 JHT=IPDE,IPDF
                      LPA=MIN(LPA,LINC(JHT))
  155             CONTINUE
  126         CONTINUE
      DO 127 JPA=IPRE,IDER
c       ... on le met dans ldeb (segment llign) ...
      LDEB(JPA-IPRE+1)=LPA
c       ... nna = longueur de la "skyline" ...
      NNA= JPA- LPA       +1
c       ... njmax = profil cumulé sur un noeud ...
      NJMAX=NJMAX+NNA
 127  continue
c
c    ... njtot = profil total ...
              NJTOT=NJTOT+NJMAX
c    ... ijmax = (profil / noeud) maxi ...
              IF(IJMAX.LT.NJMAX) IJMAX=NJMAX
              SEGDES,LLIGN
              IPRE=IDER+1
  100     CONTINUE
c
 999  CONTINUE
c
c ... fin de la boucle sur les noeuds ...
      SEGSUP TRATRA
C
C  ****  ON REPREND TOUTE LES MATRICES CONTENANT LES MULTIPLICATEURS
C  ****  POUR MULTIPLIER TOUS LEURS TERMES PAR UNE NORME ATTACHEE
C  ****  A CHAQUE MULTIPLICATEUR.PUIS ON LES ASSEMBLE.
C
*  d'abord etablir une norme generale pour le cas ou on n'arrive pas
*  a calculer la norme particuliere
      DMAXGE=xpetit
      DO 378 I=1,INC
**    write (6,*) ' assem5 diag vmav ',diag(i),vmax(i)
      DMAXGE=MAX(DMAXGE,abs(vmax(i)))
  378 CONTINUE
      if (dmaxge.lt.xpetit/xzprec*10) dmaxge=1.d0
      if (iimpi.eq.1524)
     >  write (6,*) ' nb inconnues facteur multiplicatif general ',
     >  INC,DMAXGE
      IENMU=0
c ... attention ! ici commence une boucle cachée (exécutée avec un goto 375) ...
  375 IENMU1 = IENMU
      IENMU=0
c ... boucle sur les rigidités qui calcule le nombre de matrices de blocages ...
      DO 376 I=1,NNR
          IF(.NOT.INOMUL(I)) IENMU=IENMU+1
 376  CONTINUE
c
c ... s'il n'y en a pas on saute cette partie du code ...
      IF( IENMU.EQ.0) GO TO 3750
c
c ... mimik contient les noms des variables primales ...
      MIMIK=IIMIK
      SEGACT,MIMIK
c ... boucle sur les rigidités ...
      DO 11 I=1,NNR
c    ... si celle si n'est pas une matrice de blocage on passe à la suivante ...
          IF(INOMUL(I)) GO TO 11
c
          DESCR=IRIGEL(3,I)
c    ... n3 = nb de variables primales ...
          N3=LISINC(/2)
          COER=COERIG(I)
          MELEME=IRIGEL(1,I)
          INCTRA=INCTRR(I)
          XMATRI=IRIGEL(4,I)
c    ... n2 = nombre d'éléments dans le support géométrique ...
          N2=NUM(/2)
c    ... À quoi correspond ce cas ? (pas de matrices élémentaires) ...
          IF (RE(/3).EQ.0) THEN
              INOMUL(I)=.TRUE.
              GOTO 11
          ENDIF
*          XMATRI=IMATTT(1)
*          SEGACT,XMATRI
c    ... n1 = nombre de variables duales ...
c    ... pourquoi va-t-on chercher n3 dans descr et n1 dans re ? ...
          N1=RE(/1)
          SEGINI,RA
c    ... boucle sur les éléments ...
          DO 14 IEL=1,N2
c       ... boucle sur les inconnues ...
              DO 15 ICO=1,N3
c          ... ija = n° local du noeud-support de l'inconnue n° ico ...
                  IJA=INUINV(NUM(NOELEP(ICO),IEL))
c          ... ijb = n° de l'inconnue ...
                  IJB=INCTRA(ICO)
c          ... ival contient la correspondance entre le n° local du ddl
c              et le n° d'équation (de ligne) correspondant ...
                  IVAL(ICO)=INCPO(IJB,IJA)
   15         CONTINUE
 
               DDDD= xpetit
c       ... dmax = max des valeurs absolues des termes diagonaux correspondants
c           au ddl de l'élément n° iel ...
c           COER a deja ete applique sur vmax
              DMAX=DDDD
*  la boucle suivante demarre 3 car les deux premiers sont les multiplicateurs de lagrange
      DO 19 ICO=3,N3
      DMAX=MAX(DMAX,vmax(IVAL(ICO)))
   19 CONTINUE
**    write (6,*) ' assem2 dmax dmaxge ',dmax,dmaxge
C   AUX FINS D'EVITER DES PROBLEMES DANS LA DECOMPOSITION
      IF( IIMPI. EQ.1524 ) WRITE(IOIMP,7391)DMAX,DDDD,IENMU,IENMU1
     1,I,IEL
 7391 FORMAT(' DMAX DDDD IENMU IENMU1 I  IEL',2E12.5,4I3)
**    write (6,*) ' assem2 dmax dmxge',dmax,dmaxge
      IF(DMAX.LE.xzprec*dmaxge.AND.IENMU.NE.IENMU1.AND.IEL.EQ.1)GOTO 377
      IF(DMAX.LE.xzprec*dmaxge) DMAX = DMAXGE
*  facteur de normalisation cf PV pour ne pas avoir de pivot nul
      DMAX=DMAX*1.5
*  on penalise la matrice en cas de resolution iterative
**pv  if (nucrou.eq.1) DMAX=DMAX*1D5
**    write (6,*) ' assem2 i iel dmax ',i,iel,dmax
*              XMATRI=IMATTT(IEL)
*              SEGACT,XMATRI
c       ... dmaxy = maximum des termes dans la première colonne (les 2 premiers exclus) ...
        DMAXY=SQRT(XPETIT)*1D5
*        if (norinc.eq.0) dmaxy=1.d0
*  demarrage a 3 aussi. On a toujours 1 -1 sur les LX
              DO 821 ICO=3,N1
                  DMAXY = MAX ( DMAXY, ABS( RE(ICO,1,iel)))
821           CONTINUE
c
              IF( IIMPI. EQ.9022 ) WRITE(IOIMP,7398) DMAX
 7398         FORMAT('FACTEUR MULTIPLICATIF DE NORME=',E12.5)
** coer est dans dmax, mais pas dans dmaxy
      DMAX = DMAX / DMAXY       
 
c
c       ... copie de la matrice élémentaire fois dmax*coer dans ra ...
              DO 21 ICO=1,N1
                  DO 2110 IKO=1,N1
                      RA(ICO,IKO)=RE(ICO,IKO,iel)*coer*DMAX  
 2110             CONTINUE
 21           CONTINUE
c       ... la sous-matrice de dimension 2 est multipliée par dmaxy ...
*       commenté car pose un problème dans excite et ne sert sans doute à rien
**            write (6,*) ' dmaxy dans assem5 ', dmaxy
            if (norinc.eq.0.d0) dmaxy=dmaxy*2.d0
            RA(1,1)=RA(1,1)*DMAXY
            RA(2,1)=RA(2,1)*DMAXY
            RA(1,2)=RA(1,2)*DMAXY
            RA(2,2)=RA(2,2)*DMAXY
c       ... mise à dmax des dnor correspondant aux multiplicateurs ...
              DO 22 ICO=1,2
                  dnor(ival(ico))=dmax   
   22         CONTINUE
c
c       ... boucle sur les variables primales ...
              DO 24 ICO=1,N3
c          ... ino = n° local du noeud ...
                  INO=INUINV(NUM(NOELEP(ICO),IEL))
c          ... io = n° de ligne du ddl n° ico ...
                  IO=IVAL(ICO)
      if (ico.eq.1) io1=io
      if (ico.eq.2) io2=io
c          ... cas des noeuds non totalement maitres ...
                  IF(IO.GT.NBNNMA) THEN
                      DO 1161 JJ=1,NTTMAI(/1)
                          IF (NTTMAI(JJ).EQ.INO) GOTO 1171
 1161                 CONTINUE
                      GOTO 1181
 1171                 CONTINUE
                      INO=NNOE+JJ
 1181                 CONTINUE
                  ENDIF
c          ... llign contient les lignes liées aux noeud ino ...
                  LLIGN=ILIGN(INO)
                  SEGACT,LLIGN*MOD
c          ... on rajoute le terme diagonal au diag ...
                  DIAG(IO)=DIAG(IO)+RA(ICO,ICO)
c          ... immmm contient les n°s des lignes des ddl ...
                  DO 132 JLIJ=1,IMMMM(/1)
                      JLIJ1=JLIJ
c             ... on est censé trouver le bon n° de ligne et aller au 133 ...
                      IF( IMMMM(JLIJ).EQ.IO) GO TO 133
  132             CONTINUE
c
                  IF(IIMPI.EQ.1524) WRITE(IOIMP,7354)
 7354             FORMAT( ' PREMIERE ERREUR 5')
                  CALL ERREUR(5)
                  RETURN
c
  133             CONTINUE
c
c          ... deuxième niveau de boucle sur les variables primales ...
                  DO 26 IRO=1,N3
c             ... ia = n° de colonne du ddl n° iro ...
                      IA=IVAL(IRO)
c             ... cas symétrique !!! ...
                      IF(IA.GT.IO) GO TO 26
c
c             ... jld et jlt sont égaux au début et à la fin (dans xxva et linc)
c                 des données relatives au ddl n° iro ...
                      JLT=IPPO(JLIJ1+1)
                      JLD=IPPO(JLIJ1)+1
c
c             ... on cherche le bon n° de colonne dans linc ...
                      DO 134 JL=JLD,JLT
                          JL1=JL
                          IF(LINC(JL).EQ.IA) GO TO 135
 134                  CONTINUE
c
                      IF(IIMPI.EQ.1524) WRITE(IOIMP,7355)
 7355                 FORMAT( ' DEUXIEME ERREUR 5')
                      CALL ERREUR(5)
                      RETURN
c
 135                  CONTINUE
c             ... on rajoute le coefficient dans ra au xxva correspondant ...
                      XXVA(JL1)=XXVA(JL1)+RA(ICO,IRO)
   26             CONTINUE
                  SEGDES,LLIGN
   24         CONTINUE
*  on stocke dans ittr les couples de LX
      ittr(io1)=io2
      ittr(io2)=io1
*              SEGDES,XMATRI
   14     CONTINUE
c    ... après son assemblage, une matrice de blocage n'est plus stigmatisée comme telle ...
          INOMUL(I)=.TRUE.
  377     CONTINUE
          SEGSUP,RA
   11 CONTINUE
c ... fin de la boucle sur les rigidités ...
      GO TO 375
c ... fin de la boucle cachée ...
 3750 CONTINUE
c
      DO 18 IK=1,NNR
          INCTRA=INCTRR(IK)
          SEGSUP,INCTRA
   18 CONTINUE
*  pv on desactive tout
      NNR=IRIGEL(/2)
      DO 2 IRI=1,NNR
          DESCR=IRIGEL(3,IRI)
          SEGDES DESCR
          IPT1=IRIGEL(1,IRI)
          SEGDES IPT1
          xMATRI=IRIGEL(4,IRI)
          SEGDES xMATRI
   2  CONTINUE
      INTERR(1)=NJTOT
c
      IF(IIMPI.GE.1) WRITE(IOIMP,4821) LLVNUL,NJTOT
 4821 FORMAT('Nbre de valeurs non nulles à l assemblage      ',I9,/
     &       'Nbre de valeurs à l assemblage                 ',I9)
c
c
      SEGDES,SNTT
      SEGSUP,INCTRR
      SEGDES,MDIAG
      SEGDES,MIMIK
      SEGDES,MDNOR
      SEGSUP,ITOPO
      SEGSUP,IITOP
      SEGDES,MILIGN
      SEGSUP,INUINV
      SEGDES,MMATRI
      MMMTRI=MMATRI
      SEGDES,MINCPO
      SEGSUP,IVAL,ITRA,JNOMUL,vmax
      IF(NORINC.NE.0) THEN
          CALL ASSEM0(MRIGID,2,INWUIT)
          SEGDES,MRIGID
      ELSE
          SEGDES,MRIGID
      ENDIF
      RETURN
      END
c
c
c
c
c