Télécharger ldmt2.eso

Retour à la liste

ldmt2
C LDMT2     SOURCE    GOUNAND   24/11/12    21:15:06     12076          
      SUBROUTINE LDMT2(ITRAV1,ITOPO1,INUIN1,IMINI1,MMMTRI,IPO1,INCTR1
     $     ,INCTR2,IITOP1,TRSUP,MENAGE,LDIAG,IITOPB,ITOPOB,IPODD,njtot
     $     ,INORMU)
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  ****  SUBROUTINE POUR FAIRE L'ASSEMBLAGE DE MATRICES
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)
C  ****  INORMU : =0 si on ne veut pas normaliser les multiplicateurs
C                 de Lagrange, <>0 sinon
C
C   Appelée par : LDMT1
C
C   Auteur : Michel BULIK
C
C   Date : Printemps '95
C
C   Langage : ESOPE + FORTRAN77
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
-INC CCREEL
 
-INC PPARAM
-INC CCOPTIO
-INC SMELEME
-INC SMRIGID
-INC SMMATRI
 
      SEGMENT,INUINV(NNGLOB)
      SEGMENT,ITOPO(IENNO)
      SEGMENT,ITOPOB(IENNO)
      SEGMENT,IITOP(NNOE+1)
      SEGMENT,IITOPB(NNOE+1)
      SEGMENT,IMINI(INC)
      SEGMENT,IPOS(NNOE1)
      SEGMENT,IPOD(NNOE1)
      SEGMENT,INCTRR(NIRI)
      SEGMENT,INCTRD(NIRI)
      SEGMENT,INCTRA(NLIGRE)
      SEGMENT,INCTRB(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
 
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
 
      SEGMENT,RA(N1,N1)*D
      SEGMENT JNOMUL
         LOGICAL INOMUL(NNR)
      ENDSEGMENT
      REAL*8 DMAX,COER,DDDD,DMAXY
      LOGICAL NOMUL
      LOGICAL TRSUP,MENAGE,LDIAG
 
 
      SAVE IPASG
      DATA IPASG/0/
 
 
**    write(6,*) ' ldmt2 inormu ',inormu
C
C en cas de normalisation des variables.
C
C ... On a vérifié dans RESOU que dans le cas non-symétrique soit NORINC=NORIND=0,
C     soit les deux sont non-nuls, mais on se méfie quand même ...
       inwuit=0
C       write (6,*) ' NORINC NORIND LDIAG',norinc,norind,ldiag
      IF(NORINC.NE.0.AND.NORIND.NE.0 .AND. LDIAG) then
         CALL ASSEM0(ITRAV1,1,inwuit)
      endif
C
C  PV  ON ACTIVE UNE FOIS POUR TOUTES LES MELEME DESCR... DE LA RIGIDITE
C      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
C ... dans le cas asymétrique on se doit de faire la même chose avec les vd ...
      INCTRD=INCTR2
      SEGACT,INCTRD
 
      MRIGID=ITRAV1
      SEGACT,MRIGID
      NNR=IRIGEL(/2)
C ... NNR = équivalent à NRIGEL (nb. de rigidités élémentaires) ...
      NNN=0
      SEGINI JNOMUL
      DO 1 IRI=1,NNR
         INCTRA=INCTRR(IRI)
         SEGACT INCTRA
C ... dans le cas asymétrique on se doit de faire la même chose avec les vd ...
         INCTRB=INCTRD(IRI)
         SEGACT INCTRB
         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
         if (.not.trsup) then
C    ... NA = nb de variables primales ...
            NA=LISINC(/2)
         else
C    ... NA = nb de variables duales ...
            NA=LISDUA(/2)
         endif
C    ... NNN = max des nb de variables primales ...
         NNN=MAX(NA,NNN)
C    ... INOMUL (NOn MULtiplicateurs ?) dit s'il ne s'agit pas des CL ...
         INOMUL(IRI)=.TRUE.
         IF(IPT1.ITYPEL.EQ.49) INOMUL(IRI)=.FALSE.
   1  CONTINUE
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
      IF(.NOT.LDIAG)   THEN
         ITOPO=ITOPO1
         IITOP=IITOP1
         IPOS=IPO1
      ELSE
         ITOPO=ITOPOB
         IITOP=IITOPB
         IPOS = IPODD
      ENDIF
      SEGACT,IITOP
      SEGACT,ITOPO
      SEGACT,IPOS
C       write(6,*) 'ipos', ( ipos(IU),IU=1,ipos(/1))
C       write(6,*) 'itopo', ( itopo(IU),IU=1,itopo(/1))
C       write(6,*) 'iitop', ( iitop(IU),IU=1,iitop(/1))
      INUINV=INUIN1
      SEGACT,INUINV
      SEGACT,IPOS
C ... NNOE = nombre de noeuds impliqués ...
      NNOE=IPOS(/1)-1
C ... INC = dimension de la matrice ...
      INC=IPOS(NNOE+1)
      MMATRI=MMMTRI
      SEGACT,MMATRI*MOD
      IF(LDIAG) THEN
         SEGINI,MDIAG
         IDIAG=MDIAG
      ELSE
         MDIAG=IDIAG
         SEGACT,MDIAG
      ENDIF
      SEGINI,MILIGN
      IF(TRSUP) THEN
         IILIGS=MILIGN
      ELSE
         IILIGN=MILIGN
      ENDIF
      MINCPO=IINCPO
      SEGACT,MINCPO
C ... INCDIF = nb de variables primales différentes ...
      INCDIF=INCPO(/1)
      MIPO1=IDUAPO
      SEGACT,MIPO1
casym(OK) ... on va prendre en compte aussi les vd pour déterminer
c         la taille de TRATRA ...
      INCDID=MIPO1.INCPO(/1)
      INCDIF=MAX(INCDIF,INCDID)
*
      MIMIK=IIMIK
      SEGACT,MIMIK
C ... IPV est de taille NNOE ...
      SEGINI IPV
      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       ... IIU = pointeur dans IITOP ...
            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       ... On va chercher ce maillage dans IRIGEL ...
            meleme=IRIGEL(2,IRI)
            if (meleme.eq.0) meleme=IRIGEL(1,IRI)
C       ... Boucle sur les noeuds de cet élément ...
            DO 83 I=1,NUM(/1)
C          ... IP = numéro local du I-ème noeud de l'élément ...
               IP=INUINV(NUM(I,IEL))
C          ... La première condition = les "connections" sont symétriques ...
               IF (IP.GT.INO) GOTO 83
C          ... Le noeud IP a déjà vu ce noeud ...
               IF (IPV(IP).EQ.INO) GOTO 83
C          ... Ces deux opérations se font si (IP <= INO et IPV(IP) != INO) ...
C          ... IPV = N° max des noeuds "connectés" au noeud IP ...
               IPV(IP)=INO
C          ... ICOMPT = nombre de noeuds connectés (de numéro < ou =) au noeud INO ...
               ICOMPT=ICOMPT+1
  83        CONTINUE
 81      CONTINUE
C    ... INCRED = maximum des ICOMPT sur tous les noeuds impliqués ...
         INCRED=MAX(INCRED,ICOMPT)
cdbg      write(*,*) 'LDMT2 : Noeud ',INO,', ICOMPT = ',ICOMPT,
cdbg     &           ' -> INCRED = ',INCRED
  80  CONTINUE
cdbg      segprt,ipv
      SEGSUP IPV
*
      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     -----------------------------
 
      LLVNUL=0
C ... On n'initialise plus IJMAX ici - ceci se fait un niveau au dessus ...
C-OLD      IJMAX=0
C ... La taille de MDNOR est INC (nombre d'équations) ...
      MDNOR=0
      IF (LDIAG) THEN
         SEGINI MDNOR
         IDNORM=MDNOR
C ... On a vérifié dans RESOU que dans le cas non-symétrique soit NORINC=NORIND=0,
C     soit les deux sont non-nuls ...
         SEGINI,MDNO1
         IDNORD=MDNO1
         DO 20 I=1,INC
            DNOR(I)=1.D0
            MDNO1.DNOR(I)=1.D0
   20    CONTINUE
      else
         mdnor=idnorm
         segact mdnor
         mdno1=idnord
         segact mdno1
      ENDIF
 
C ... Si une normalisation est imposée, les valeurs dans MDNOR et MDNO1 seront modifiées par ASNS3 ...
      midua=iidua
      IF(NORINC.NE.0 .AND. NORIND.NE.0 .AND. LDIAG)
     &             CALL ASNS3(MDNOR,MIMIK,MINCPO,mdno1,midua,mipo1)
 
C  ****  BOUCLE   *100* SUR LES NUMEROS DE NOEUDS QUE L'ON ASSEMBLE
C        -------------------------------------------------------------
 
      DO 100 INO=1,NNOE
cdbg      write(*,*) 'LDMT2: ----------------------------------'
cdbg      write(*,*) 'LDMT2: Noeud N° ',INO
         DO 101 IIT=1,INCDIF
            MTRA(IIT)=0
 101     CONTINUE
 
C ... Attention à l'égalité des IPOS et IPOD qui est testée dans ASNS1 ...
 
C    ... IPRE = N° de la première ligne (ou colonne) concernée par le noeud INO ...
         IPRE=IPOS(INO  )+1
C    ... IDER = N° de la dernière ligne (ou colonne) concernée par le noeud INO ...
         IDER=IPOS(INO+1)
cdbg      write(*,*) 'LDMT2:    Lignes de ',IPRE,' jusqu''à ',IDER
         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
C    ... MAXELE = nb d'éléments qui contiennent le noeud No INO ...
         MAXELE= (IITOP(INO+1) -IITOP(INO))/2
         DO 99 IELE=1,MAXELE
cdbg      write(*,*) 'LDMT2:    - - - - - - - - - - - -'
cdbg      write(*,*) 'LDMT2:    Elément N° ',IELE
            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 ...
casym(OK) ... ici il faudra choisir en fonction de TRSUP si prendre INCTRR ou INCTRD ...
            INCTRA=INCTRR(IRI)
            INCTRB=INCTRD(IRI)
C       ... xMATRI = tableau de pointeurs vers les matrices élémentaires ...
            xMATRI=IRIGEL(4,IRI)
C       ... COER = coefficient multiplicateur ...
            COER=COERIG(IRI)
C       ... NOn MULtiplicateur ? ...
            NOMUL=INOMUL(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 ... Maintenant, dans la boucle 96 (nouvelle par rapport à ASSEM2) on
C     va parcourir les lignes (colonnes) pour trouver celles, qui
C     s'appuyent sur le noeud INO, cas inférieur (supérieur) ...
 
            NA=0
            IF(TRSUP) THEN
               NIN=LISINC(/2)
            ELSE
               NIN=LISDUA(/2)
            ENDIF
            DO 96 ICO=1,NIN
               IF(TRSUP) THEN
                  IJA=INUINV(NUM(NOELEP(ICO),IEL))
                  IJB=INCTRA(ICO)
               ELSE
                  IJA=INUINV(NUM(NOELED(ICO),IEL))
                  IJB=INCTRB(ICO)
               ENDIF
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
cdbg      write(*,*) 'LDMT2:         IJA == INO => NA devient ',NA
                  ITRA(NA,1)=ICO
                  IF(TRSUP) THEN
                     ITRA(NA,2)=INCPO(IJB,IJA)
                  ELSE
                     ITRA(NA,2)=MIPO1.INCPO(IJB,IJA)
                  ENDIF
               ENDIF
 96         CONTINUE
 
C ... Dans la boucle 98 on va remplir le tableau IVAL avec les numéros
C     de colonnes (lignes) pour tous les DDLs primaux (duaux) de la matrice,
C     ceci concerne le triangle inférieur (supérieur) ...
casym(OK) ... la boucle devra se faire soit sur les duales soit sur les primales ...
 
            IF(TRSUP) THEN
               NIN=LISDUA(/2)
            ELSE
               NIN=LISINC(/2)
            ENDIF
C       ... Boucle sur les variables primales (duales) de l'élément ...
            DO 98 ICO=1,NIN
               IF(TRSUP) THEN
                  IJA=INUINV(NUM(NOELED(ICO),IEL))
                  IJB=INCTRB(ICO)
                  IVAL(ICO)=MIPO1.INCPO(IJB,IJA)
cdbg      write(*,*) 'LDMT2:      vd locale N°',ICO,'-> vd globale N°',IJB
cdbg      write(*,*) 'LDMT2:      => N° d''équation = ',IVAL(ICO)
               ELSE
C          ... IJA = numéro local du noeud-support du DDL ...
casym(OK) ... ici il faudra choisir en fonction de TRSUP si prendre NOELEP ou NOELED ...
                  IJA=INUINV(NUM(NOELEP(ICO),IEL))
C             ... IJB = indice du DDL ...
casym(OK) ... INCTRA ou INCTRB ??? ...
                  IJB=INCTRA(ICO)
C          ... On met dans IVAL(ICO) le N° de la colonne correspondant au noeud IJA
C              et DDL No IJB ...
casym(OK) ... ici il faudra choisir en fonction de TRSUP si prendre MINCPO ou MIPO1 (?) ...
                  IVAL(ICO)=INCPO(IJB,IJA)
cdbg      write(*,*) 'LDMT2:      vp locale N°',ICO,'-> vp globale N°',IJB
cdbg      write(*,*) 'LDMT2:      => N° de la colonne = ',IVAL(ICO)
               ENDIF
  98        CONTINUE
 
C            XMATRI=IMATTT(IEL)
C            SEGACT,XMATRI
cdbg      segprt,xmatri
C
C  ****  BOUCLE  *95* SUR LES INCONNUES DE LA PETITE MATRICE
C
            DO 95 INCC=1,NA
C          ... INCO = le N° de la ligne (cas inf.) ou colonne (cas sup.)
C              du DDL dual (primal) N° INCC du noeud INO ...
casym(OK) ... ici il faudra choisir en fonction de TRSUP quoi prendre ...
               INCO=ITRA(INCC,2)
C          ... Boucle sur les DDL's du noeud INO ...
               DO 90 IK=1,NIN
C             ... IO = N° de la colonne (la ligne) correspondant au DDL N° IK de l'élément
casym(OK) ... ici il faudra choisir en fonction de TRSUP quoi prendre ...
                  IO=IVAL(IK)
C             ... Si IO > INCO on ne fait plus rien dans cette boucle ...
C             ... Car on ne parcourt que la moitié de la matrice !!! ...
                  IF(IO.GT.INCO)  GO TO 90
C             ... ILOC = N° de ligne relatif / IPRE du noeud INO ...
                  ILOC=INCO-IPRE+1
C             ... JJ = N° (dans l'élément) de la ligne (de la colonne) correspondant au
C                 DDL N° INCC du noeud INO ...
casym(OK) ... attention à ce qui se passe ici ... !!!
                  JJ=ITRA(INCC,1)
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 ...
C             ... ILTT doit contenir l'indice IMMTT en fonction du N° de la colonne; s'il est
C                 nul, le terme n'a pas encore été assemblé, il faudra donc initialiser ...
                  ILTT= LTRA(IO,ILOC)
                  IF(ILTT.EQ.0) THEN
C                ... LLVVA est initialisé à 0 au début de la boucle 100 et sert à compter les
C                    termes assemblés ...
                     LLVVA=LLVVA+1
C                ... Les MTRA sont mis à 0 au début de la boucle 100 ...
C                ... Ici on incrémente MTRA(ILOC) et on met la nouvelle valeur dans IMMTT,
C                    qui sera le premier indice dans différents tableaux de TRATRA ...
                     IMMTT=MTRA(ILOC)+1
                     MTRA(ILOC)=IMMTT
C                ... puis on initialise le nouveau terme ...
                     XTRA(IMMTT,ILOC)=0.D0
C                ... son N° de colonne (cas inf.) ou ligne (cas supérieur) ...
                     NTRA(IMMTT,ILOC)=IO
C                ... et l'information inverse ...
                     LTRA(IO,ILOC)=IMMTT
                     ILTT=IMMTT
                  ENDIF
C             ... Si ce n'est pas une condition aux limites, on assemble ...
                  IF(NOMUL) THEN
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 ...
                     IF(TRSUP) THEN
C                write(6,*) ' iltt iloc ' , iltt,iloc
C           write(6,*) 'ik jj iel re coer ',ik,jj,iel,re(ik,jj,iel),coer
                        XTRA(ILTT,ILOC)=XTRA(ILTT,ILOC)+RE(IK,JJ,iel)
     $                       *COER
cdbg      write(*,*) 'LDMT2: On rajoute le terme L',IK,', C',JJ,' soit',
cdbg     &           RE(IK,JJ),' à la col',IPRE+ILOC-1,' et lig',
cdbg     &           NTRA(ILTT,ILOC)
                     ELSE
C                     write(6,*) ' iltt iloc ' , iltt,iloc
C        write(6,*) 'ik jj iel re coer ',ik,jj,iel,re(jj,ik,iel),coer
                        XTRA(ILTT,ILOC)=XTRA(ILTT,ILOC)+RE(JJ,IK,iel)
     $                       *COER
cdbg      write(*,*) 'LDMT2: On rajoute le terme L',JJ,', C',IK,' soit',
cdbg     &           RE(JJ,IK),' à la lig',IPRE+ILOC-1,' et col',
cdbg     &           NTRA(ILTT,ILOC)
                     ENDIF
                  ENDIF
   90          CONTINUE
   95       CONTINUE
C            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
         SEGINI,LLIGN
         ILIGN(INO)=LLIGN
         NBA=0
C         write(6,*) ' ider ipre llvnul ino ',ider,ipre,llvnul,ino
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       ... que l'on met au bon endroit dans IMMMM (faisant partie de LLIGN) ...
            IMMMM(JPA)=IIIN
C       ... IPPO(I)+1 = adresse du début (dans XXVA et LINC) des données relatives au DDL N° I ...
            IPPO(JPA)=NBA
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               write(6,*) ' nba ipak jpa, ', nba ,ipak,jpa
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(LDIAG) THEN
C               write(6,*) ' iiin  nba xxva(nba)', iiin,nba,xxva(nba)
                  IF(IIIN.EQ.IUNPAK) DIAG(IIIN)=XXVA(NBA)
               ENDIF
  121       CONTINUE
  120    CONTINUE
C    ... Le pointeur vers la fin de la dernière ligne ...
         IPPO(NA+1)= NBA
         NJMAX= 0
C  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)=INO
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    ... NJTOT = profil total ...
         NJTOT=NJTOT+NJMAX
C    ... IJMAX = (profil / noeud) maxi ...
         IF(IJMAX.LT.NJMAX) IJMAX=NJMAX
         SEGDES,LLIGN
  100 CONTINUE
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,*) ' assem2 diag vmav ',diag(i),vmax(i)
**    vmax(i)=max(vmax(i),abs(diag(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*10) dmaxge=1.d0
      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 ... S'il n'y en a pas on saute cette partie du code ...
      IF( IENMU.EQ.0) GO TO 3750
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
         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
C         XMATRI=IMATTT(1)
C         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
C           JCARDO => pour les fous qui veulent se passer de la
C                     normalisation des mult. de Lagrange :)
            IF (INORMU.EQ.0) THEN
               DMAX=1.D0
               DMAXY=DMAX
            ELSE
               DDDD= xpetit
C           ... DMAX = max des valeurs absolues des termes diagonaux correspondants
C               au DDL de l'élément N° IEL ...
               DMAX=xpetit
               IF((.NOT.LDIAG).AND.(IDIAG.EQ.0)) THEN
C          write(*,*) 'Erreur interne dans LDMT2 : MDIAG pas encore',
C         &           ' initialisé !!!'
                  RETURN
               ENDIF
*  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)),abs(diag(ival(ico))))
 19            CONTINUE
 
C           ... AUX FINS D'EVITER DES PROBLEMES DANS LA DECOMPOSITION
               if(iimpi.eq.1524) WRITE(IOIMP,7391) DMAX,DDDD,IENMU,
     &                                            IENMU1,I,IEL
 7391          FORMAT(' DMAX DDDD IENMU IENMU1 I  IEL',2E12.5,4I3)
 
               IF(DMAX.LE.xpetit.AND.IENMU.NE.IENMU1
     &             .AND.IEL.EQ.1) GOTO 377
               IF(DMAX.LE.xzprec*dmaxge) DMAX = DMAXGE
               DMAX=DMAX*1.5D0
C                XMATRI=IMATTT(IEL)
C                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 ICO=3,N1
                  DMAXY = MAX ( DMAXY, ABS( RE(ICO,1,iel)))
               enddo
*  demarrage a 3 aussi. On a toujours 1 -1 sur les LX
               DO iko=3,N3
                  DMAXY = MAX ( DMAXY, ABS( RE(1,iko,iel)))
               enddo
               DMAX = DMAX / DMAXY
*    Attention
*    on reprend le dmax du premier passage pour appliquer le meme facteur sur les lignes et les colonnes
*    car le max de la ligne n'est pas forcement celui de la colinne
               if (.not.ldiag) dmax=dnor(ival(1))
            ENDIF
            IF( IIMPI. EQ.1524 ) WRITE(IOIMP,7398) DMAX
 7398       FORMAT('  facteur multiplicatif de norme ',e12.5)
C       ... Copie de la matrice élémentaire fois DMAX*COER dans RA ...
            DO 21 ICO=1,N1
               DO 2110 IKO=1,N3
                  RA(ICO,IKO)=RE(ICO,IKO,iel)*DMAX*COER
 2110          CONTINUE
 21         CONTINUE
C       ... La sous-matrice de dimension 2 est multipliée par DMAXY ...
**  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
C       ... Mise à DMAX des DNOR correspondant aux multiplicateurs ...
            IF(LDIAG) THEN
               DO 22 ICO=1,2
                  DNOR(IVAL(ICO))=DMAX
**                write (6,*) 'il dnor ',ival(ico),dnor(ival(ico))
                  mdno1.DNOR(IVAL(ICO))=DMAX
   22          CONTINUE
            ENDIF
CMB     ... Je suppose que l'on s'en servira pour multiplier les composantes
C           du second membre pour que les déplacements imposés soient corrects ...
C       ... Boucle sur les variables primales ...
**             write(6,*) 'ldmt2 trsup ',trsup
            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          ... LLIGN contient les lignes liées aux noeud INO ...
               LLIGN=ILIGN(INO)
               SEGACT,LLIGN*MOD
C          ... On rajoute le terme diagonal au DIAG ...
               IF(LDIAG) 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
               IF(IIMPI.EQ.1524) WRITE(IOIMP,7354)
 7354          FORMAT( ' PREMIERE ERREUR 5')
               CALL ERREUR(5)
               RETURN
  133          CONTINUE
C          ... Deuxième niveau de boucle sur les variables primales ...
               DO 26 IRO=1,N1
C             ... IA = N° de colonne du DDL N° IRO ...
                  IA=IVAL(IRO)
C             ... Cas symétrique !!! ...
                  IF(IA.GT.IO) GO TO 26
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             ... 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
                  IF(IIMPI.EQ.1524) WRITE(IOIMP,7355)
 7355             FORMAT( ' DEUXIEME ERREUR 5')
                  CALL ERREUR(5)
                  RETURN
 135              CONTINUE
C             ... On rajoute le coefficient dans RA au XXVA correspondant ...
                  if(.not.trsup) then
                     XXVA(JL1)=XXVA(JL1)+RA(ICO,IRO)
                  else
                     XXVA(JL1)=XXVA(JL1)+RA(IRO,ICO)
                  endif
   26          CONTINUE
               SEGDES,LLIGN
   24       CONTINUE
C  on stocke dans ittr les couples de LX
            ittr(io1)=io2
            ittr(io2)=io1
C            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
**    do i = 1,inc
**     if (ittr(i).eq.0.and.diag(i).ne.vmax(i))
**   >   write (6,*) 'diag vmax ',i,diag(i),vmax(i)
**    enddo
C ... C'est l'ménaage finaale ...
      DO 18 IK=1,NNR
         INCTRA=INCTRR(IK)
         IF(MENAGE) THEN
            SEGSUP,INCTRA
         ELSE
            SEGDES,INCTRA
         ENDIF
 18   CONTINUE
C  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-old      IF(IPASG.EQ.0) CALL ERREUR(-278)
      IF(MENAGE.AND.IPASG.EQ.0) CALL ERREUR(-278)
      IPASG=1
 
      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)
         SEGSUP,IITOP,ITOPO
      IF(MENAGE) THEN
C     SEGSUP,IITOP
ccc   SEGSUP,IMINI
         SEGSUP,INUINV
C         SEGSUP,ITOPO
         SEGSUP,INCTRR,INCTRD
      ELSE
C         SEGDES,IITOP
ccc         SEGDES,IMINI
c         SEGDES,INUINV
C         SEGDES,ITOPO
c         SEGDES,INCTRR
      ENDIF
      IF (MENAGE) THEN
         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
         segact mrigid*nomod
      ENDIF
C      write(6,*) ' diag'
C      write(6,*) ( diag(iou),iou=1,diag(/1))
      SEGDES,MDIAG
      SEGDES,MIMIK
      IF (MDNOR.NE.0) SEGDES,MDNOR
      SEGDES,MILIGN
      SEGDES,MMATRI
      MMMTRI=MMATRI
      SEGDES,MINCPO
      SEGDES,MIPO1
      SEGDES,IPOS
      SEGSUP,IVAL,ITRA,JNOMUL,vmax
      IF(NORINC.NE.0 .AND. MENAGE) THEN
         CALL ASSEM0(MRIGID,2,inwuit)
         SEGDES,MRIGID
      ELSE
         SEGDES,MRIGID
      ENDIF
      RETURN
      END
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales