C ASSEM2    SOURCE    PV090527  26/04/30    21:15:08     12529          
      SUBROUTINE ASSEM2(ITRAV1,ITOPO1,INUIN1,IMINI1,MMMTRI,IPO1,INCTR1
     $     ,IITOP1)
C
C  ****  SUBROUTINE POUR FAIRE L'ASSEMBLAGE DE MATRICES SYMETRIQUES
C        EN VUE D'UN TRAITEMENT PAR METHODE DE KROUT.
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  ****  IMINI1 : POINTEUR SEGMENT DE TRAVAIL IMINI (VOIR ASSEM1)
C  ****  IPO1   : POINTEUR SEGMENT DE TRAVAIL IPOS  (VOIR ASSEM1)
C  ****  MMMTRI : POINTEUR OBJET MATRICE TRIANGULARISEE (NON MODIFIE)
C                 (VOIR SMMATRI)
*
*  modif janvier 2015 toutes les inconnues d'un noeud commencent à la même colonne
*  modif mars 2018 la verification de symetrie des matrices elementaires est externalisee
*
      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,IMINI(INC)
      SEGMENT,IPOS(NNOE1)
      SEGMENT,INCTRR(NIRI)
      SEGMENT,INCTRA(NLIGRE)
      SEGMENT,IPV(NNOE)
      SEGMENT,VMAX(INC)

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  ****  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
      SEGMENT,RA(N1,N1)*D
      SEGMENT JNOMUL
         LOGICAL INOMUL(NNR)
      ENDSEGMENT
      REAL*8 DMAX,COER,DDDD,DMAXY,DMAXGE
      LOGICAL NOMUL
**    nbver=0
*
*  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
      INCTRR=INCTR1
      SEGACT,INCTRR

      MRIGID=ITRAV1
      SEGACT,MRIGID
      NNR=IRIGEL(/2)
      NNN=0
      SEGINI JNOMUL
      DO 1 IRI=1,NNR
         INCTRA=INCTRR(IRI)
         SEGACT INCTRA
         DESCR=IRIGEL(3,IRI)
         SEGACT, DESCR
         IPT1=IRIGEL(1,IRI)
         SEGACT IPT1
         ipt2=IRIGEL(2,IRI)
         if (ipt2.ne.0) segact,ipt2
*      XMATRI=IRIGEL(4,IRI)
*      SEGACT XMATRI
         NA=LISINC(/2)
         NNN=MAX(NA,NNN)
         INOMUL(IRI)=.TRUE.
         IF(IPT1.ITYPEL.EQ.49) INOMUL(IRI)=.FALSE.
   1  CONTINUE
      SEGINI,IVAL
      SEGINI,ITRA
C
C  ****  ACTIVATION DES SEGMENTS DE TRAVAILS ET DE MMATRI
C
*      IMINI=IMINI1
*      SEGACT,IMINI
*      N1=IMINI(/1)
      ITOPO=ITOPO1
      SEGACT,ITOPO
      IITOP=IITOP1
      SEGACT,IITOP
*
      INUINV=INUIN1
      SEGACT,INUINV
      IPOS=IPO1
      SEGACT,IPOS
*
      NNOE=IPOS(/1)-1
*
      INC=IPOS(NNOE+1)
      MMATRI=MMMTRI
      SEGACT,MMATRI*MOD
      SEGINI,MDIAG
      SEGINI,MILIGN
      IILIGN=MILIGN
      IDIAG=MDIAG
      MINCPO=IINCPO
      SEGACT,MINCPO
      INCDIF=INCPO(/1)
*
      MIMIK=IIMIK
      SEGACT,MIMIK
      SEGINI IPV
      INCRED=0
      DO 80 INO =1,NNOE
         ICOMPT=0
         MAXELE = (IITOP(INO+1)-IITOP(INO))/2
         DO 81 IELE=1,MAXELE
            IIU=IITOP(INO) + IELE + IELE -2
            IEL=ITOPO(IIU)
            IRI=ITOPO(IIU+1)
            meleme=IRIGEL(2,IRI)
            if (meleme.eq.0) meleme=IRIGEL(1,IRI)
            DO 83 I=1,NUM(/1)
               IP=INUINV(NUM(I,IEL))
               IF (IP.GT.INO) GOTO 83
               IF (IPV(IP).EQ.INO) GOTO 83
               IPV(IP)=INO
               ICOMPT=ICOMPT+1
 83         CONTINUE
 81      CONTINUE
         INCRED=MAX(INCRED,ICOMPT)
  80  CONTINUE
      SEGSUP IPV
*
      INCRED=INCRED*INCDIF
      SEGINI TRATRA
      segini vmax
C
C  ****  BOUCLE   *100* SUR LES NUMEROS DE NOEUDS QUE L'ON ASSEMBLE
C
      LLVNUL=0
      IJMAX=0
      NJTOT=0
      NNOE=IPOS(/1)-1
      SEGINI MDNOR
      IDNORM=MDNOR
*
* en cas de normalisation des variables.
*
      IF(NORINC.NE.0)  THEN
         inwuit=0
         CALL ASSE10(ITRAV1,1,MDNOR,MIMIK,MINCPO,INUIN1,inwuit)
      ELSE
         DO 53 IU=1,INC
            DNOR(IU)=1.d0
 53      CONTINUE
      ENDIF
*  verif de la symetrie des matrices elementaires (sauf le super element qui peut etre tres grand)
      do jr=1,irigel(/2)
       IPT6=IRIGEL(1,JR)
       if (ipt6.itypel.ne.28) then
        XMATRI=IRIGEL(4,jr)
        SEGACT XMATRI*mod
**      nbver=max(nbver,re(/1))
        if(symver.ne.1) call versym(re,re(/1),re(/2),re(/3),0)
        if (ierr.ne.0) return
        symre=0
        symver=1
        segdes xmatri
       endif
      enddo
*     On balaye les noeuds dans l'ordre des elements
*     (Pourquoi ? Dans le non-symétrique, ce n'est pas fait ldmt2.eso )
      DO 100 JR=1,IRIGEL(/2)
         ipt6=IRIGEL(2,JR)
         if (ipt6.eq.0) ipt6=IRIGEL(1,JR)
         DO 101 JL=1,IPT6.num(/2)
            DO 102 JP=1,IPT6.num(/1)
               INO=INUINV(IPT6.NUM(JP,JL))
**     DO 100 INO=1,NNOE
               IF (ILIGN(INO).NE.0) GOTO 102
               DO 103 IIT=1,INCDIF
                  MTRA(IIT)=0
 103           CONTINUE
               IPRE=IPOS(INO)+1
               IDER=IPOS(INO+1)
               LLVVA=0
C
C  ****  BOUCLE   *99*   SUR LES ELEMENTS TOUCHANT LE NOEUD INO
C                        POUR LES ELEMNTS MULTIPLICATEUR ON NE FAIT PAS
C                        L'ASSEMBLAGE
C
               MAXELE= (IITOP(INO+1) -IITOP(INO))/2
               DO 99 IELE=1,MAXELE
                  IIU=IITOP(INO) + IELE + IELE - 2
                  IEL=ITOPO(IIU)
                  IRI=ITOPO(IIU+1)
                  MELEME=IRIGEL(1,IRI)
                  DESCR=IRIGEL(3,IRI)
                  INCTRA=INCTRR(IRI)
                  XMATRI=IRIGEL(4,IRI)
                  SEGACT XMATRI
                  COER=COERIG(IRI)
C
C  ****  NOMUL =.FALSE.  IL EXISTE  UN MULTUIPLICATEUR
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
                  NIN=LISINC(/2)
                  NOMUL=INOMUL(IRI)
                  NA=0
                  DO 98 ICO=1,NIN
                     IJA=INUINV(NUM(NOELEP(ICO),IEL))
                     IJB=INCTRA(ICO)
                     IVAL(ICO)=INCPO(IJB,IJA)
                     IF(IJA.NE.INO)  GO TO 98
                     NA=NA+1
                     ITRA(NA,1)=ICO
                     ITRA(NA,2)=IVAL(ICO)
 98               CONTINUE
*      XMATRI=IMATTT(IEL)
*      SEGACT,XMATRI
C
C  ****  BOUCLE  *95* SUR LES INCONNUES DE LA PETITE MATRICE
C
                  DO 90 INCC=1,NA
                     INCO=ITRA(INCC,2)
                     if (inco.gt.ider) goto 90
                     ILOC=INCO-IPRE+1
                     JJ=ITRA(INCC,1)
                     DO 95 IK=1,NIN
                        IO=IVAL(IK)
                        IF(IO.GT.INCO)  GO TO 95
*?                        IPOO=IK*NIN - NIN
*?                        IPO=IPOO+JJ
                        ILTT= LTRA(IO,ILOC)
                        IF(ILTT.EQ.0) THEN
                           LLVVA=LLVVA+1
                           IMMTT=MTRA(ILOC)+1
                           MTRA(ILOC)=IMMTT
                           XTRA(IMMTT,ILOC)=0.D0
                           NTRA(IMMTT,ILOC)=IO
                           LTRA(IO,ILOC)=IMMTT
                           ILTT=IMMTT
                        ENDIF
                        IF(NOMUL) THEN
**    XTRA(ILTT,ILOC)=XTRA(ILTT,ILOC)+RE(JJ,IK,IEL)*COER
** on utilise la symetrie de re
                           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 MOUVENENT 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
               NA = IDER-IPRE+1
               LLVNUL=LLVNUL+LLVVA
               SEGINI,LLIGN
               ILIGN(INO)=LLIGN
               NBA=0
               DO 120 JPA=1,NA
                  IIIN=IPRE+JPA -1
                  IMMMM(JPA)=IIIN
                  IPPO(JPA)=NBA
                  DO 121 IPAK = 1,MTRA(JPA)
                     IUNPAK=NTRA(IPAK,JPA)
                     LTRA(IUNPAK,JPA)=0
** pas mur pour le test suivant
***   if (abs(xtra(ipak,jpa)).gt.xpetit) then
                     NBA=NBA+1
                     LINC(NBA)=IUNPAK
                     XXVA(NBA)=XTRA(IPAK,JPA)
**    write (6,*) 'assem2 iiin nba xxva(nba)',
**   >                    iiin,nba,xxva(nba)
                     vmax(iiin)=max(abs(xxva(nba)),vmax(iiin))
***   endif
                     IF(IIIN.EQ.IUNPAK) DIAG(IIIN)=xtra(ipak,jpa)
 121              CONTINUE
 120           CONTINUE
               IPPO(NA+1)= NBA
               NJMAX= 0
*  recherche du mini globale sur toutes les inconnues
               LPA=IPRE
               DO 126 JPA=IPRE,IDER
                  IPNO(JPA)=INO
                  IPDE=IPPO(JPA-IPRE+1)+1
                  IPDF=IPPO(JPA-IPRE+2)
                  DO 155 JHT=IPDE,IPDF
                     LPA=MIN(LPA,LINC(JHT))
 155              CONTINUE
 126           CONTINUE
               DO 127 JPA=IPRE,IDER
                  LDEB(JPA-IPRE+1)=LPA
                  NNA= JPA- LPA       +1
                  NJMAX=NJMAX+NNA
 127           continue
               NJTOT=NJTOT+NJMAX
               IF(IJMAX.LT.NJMAX) IJMAX=NJMAX
               SEGDES,LLIGN
 102        CONTINUE
 101     CONTINUE
 100  CONTINUE
      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,*) ' assem2 diag vmav ',diag(i),vmax(i)
         DMAXGE=MAX(DMAXGE,abs(vmax(i)))
  378 CONTINUE
      if (iimpi.ne.0   )
     >  write (6,*) ' nb inconnues facteur multiplicatif general ',
     >  INC,DMAXGE
      if (dmaxge.lt.xpetit/xzprec) dmaxge=1.d0
      IENMU=0
 375  IENMU1 = IENMU
      IENMU=0
      DO 376 I=1,NNR
         IF(.NOT.INOMUL(I)) IENMU=IENMU+1
 376  CONTINUE
      IF( IENMU.EQ.0) GO TO 3750
      MIMIK=IIMIK
      SEGACT,MIMIK
      DO 11 I=1,NNR
         IF(INOMUL(I)) GO TO 11
         DESCR=IRIGEL(3,I)
         N3=LISINC(/2)
         COER=COERIG(I)
         MELEME=IRIGEL(1,I)
         INCTRA=INCTRR(I)
         XMATRI=IRIGEL(4,I)
         SEGACT XMATRI
         N2=NUM(/2)
         IF (RE(/3).EQ.0) THEN
            INOMUL(I)=.TRUE.
            SEGDES XMATRI
            GOTO 11
         ENDIF
*      XMATRI=IMATTT(1)
*      SEGACT,XMATRI
         N1=RE(/1)
c     ERREUR 756 : La matrice de rigidite n'est pas carree
         if (n1.ne.re(/2)) call erreur(756)
         if (ierr.ne.0) return
         SEGINI,RA
         DO 14 IEL=1,N2
            DO 15 ICO=1,N3
               IJA=INUINV(NUM(NOELEP(ICO),IEL))
               IJB=INCTRA(ICO)
               IVAL(ICO)=INCPO(IJB,IJA)
 15         CONTINUE
            DMAX=xpetit
*  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,IENMU,IENMU1,I
     $           ,IEL
 7391       FORMAT(' DMAX  IENMU IENMU1 I  IEL',1E12.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.5D0
*  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
            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
*     if (dmaxy.lt.1d-50) write (6,*) (re(ico,1),ico=1,n1)
**    if (dmaxy.lt.1d+50)write (6,*)' assem2 dmax dmaxy ',
**   >       dmax,dmaxy,dmaxge
            DMAX = DMAX / DMAXY
            IF( IIMPI. EQ.1524 ) WRITE(IOIMP,7398) DMAX
 7398       FORMAT('  facteur multiplicatif de norme ',e12.5)
            DO 21 ICO=1,N1
               DO 2110 IKO=1,N1
                  RA(ICO,IKO)=RE(ICO,IKO,IEL)*coer*DMAX
 2110          CONTINUE
 21         CONTINUE
**    write (6,*) ' dmax ',dmax
**  si on ne booste pas l'egalite des mults on a des problemes de precision sur ceux ci
            if (norinc.eq.0) 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
            DO 22 ICO=1,2
               DNOR(IVAL(ICO))=DMAX
 22         CONTINUE
*      if (abs(dmax).gt. 1d50) write (6,*) ' assem2 dmax dmaxy ',dmax,
*     >  dmaxy,dmaxge
            DO 24 ICO=1,N3
               INO=INUINV(NUM(NOELEP(ICO),IEL))
               IO=IVAL(ICO)
               if (ico.eq.1) io1=io
               if (ico.eq.2) io2=io
               LLIGN=ILIGN(INO)
               SEGACT,LLIGN*MOD
               DIAG(IO)=DIAG(IO)+RA(ICO,ICO)
               DO 132 JLIJ=1,IMMMM(/1)
                  JLIJ1=JLIJ
                  IF( IMMMM(JLIJ).EQ.IO) GO TO 133
 132           CONTINUE
               IF(IIMPI.EQ.1524) WRITE(IOIMP,7354)
 7354          FORMAT( ' PREMIERE ERREUR 5')
               CALL ERREUR(5)
               RETURN
 133           CONTINUE
*      IREE=ICO*N3-ICO
               DO 26 IRO=1,N3
                  IA=IVAL(IRO)
                  IF(IA.GT.IO) GO TO 26
*      IF(IRO.LE.ICO) IRE=(ICO*(ICO-1))/2+IRO
*      IF(IRO.GT.ICO) IRE=(IRO*(IRO-1))/2+ICO
                  JLT=IPPO(JLIJ1+1)
                  JLD=IPPO(JLIJ1)+1
                  DO 134 JL=JLD,JLT
                     JL1=JL
                     IF(LINC(JL).EQ.IA) GO TO 135
 134              CONTINUE
                  IF(IIMPI.NE.1524) WRITE(IOIMP,7355)
 7355             FORMAT( ' DEUXIEME ERREUR 5')
                  CALL ERREUR(5)
                  RETURN
 135              CONTINUE
                  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
 14      CONTINUE
         INOMUL(I)=.TRUE.
 377     CONTINUE
         SEGSUP,RA
 11   CONTINUE
      GO TO 375
 3750 CONTINUE
      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)
         XMATRI=IRIGEL(4,IRI)
         SEGDES XMATRI
   2  CONTINUE

      INTERR(1)=NJTOT
      IF(IIMPI.EQ.1457) WRITE(IOIMP,4821) LLVNUL,NJTOT
 4821 FORMAT('  NB DE VALEURS NON NULLES DANS LA MATRICE ',I9,/
     #       '  NB DE VALEURS DANS LA MATRICE            ',I9)
      IF(NORINC.NE.0) THEN
         CALL ASSE10(MRIGID,2,MDNOR,MIMIK,MINCPO,INUIN1,inwuit)
          SEGDES,MRIGID
      ELSE
         SEGDES,MRIGID
      ENDIF
      SEGSUP,INCTRR
      segact mrigid*mod
      NNR=IRIGEL(/2)
      DO  IRI=1,NNR
         IPT2=IRIGEL(2,IRI)
         IF (IPT2.NE.0) THEN
            SEGSUP IPT2
            IRIGEL(2,IRI)=0
         ENDIF
      ENDDO
      segdes mrigid
      SEGDES,MDIAG
      SEGDES,MIMIK
      SEGDES,MDNOR
ccc      SEGSUP,IMINI
      SEGSUP,ITOPO
      SEGSUP,IITOP
      SEGDES,MILIGN
      SEGSUP,INUINV
      SEGDES,MMATRI
      MMMTRI=MMATRI
      SEGDES,MINCPO
      SEGDES,IPOS
      SEGSUP,IVAL,ITRA,JNOMUL,vmax
**    write(6,*) 'taille maxi matrice elementaire ',nbver
      RETURN
      END