Sources Esope | Cast3M

extr13
C EXTR13    SOURCE    PV        17/09/29    21:15:10     9578           
      SUBROUTINE EXTR13(IRIG,IRET)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
************************************************************************
*
*                             E X T R 1 3
*                             -----------
*
* FONCTION:
* ---------
*
*     Sous-programme pour extraire d'une rigidit{ :
*       - les matrices élémentaires symétriques si IRET = 3
*       - les matrices élémentaires antisymétriques si IRET = 4
*
* MODULES UTILISES:
* -----------------
*
 
-INC PPARAM
-INC CCOPTIO
-INC SMRIGID
*
* PARAMETRES:   (e)=ENTREE   (s)=SORTIE   (+ = CONTENU DANS UN COMMUN)
* -----------
*
*     IRIG    (e)  pointeur sur l'objet RIGIDITE.
*     IRET    (e)  = 3 , matrice symétrique.
*                  = 4 , matrice antisymétrique.
*
* AUTEUR, DATE DE CREATION:
* -------------------------
*
*     Lionel VIVAN       juin 1991
*
* LANGAGE:
* --------
*
*     ESOPE + FORTRAN77
*
************************************************************************
*
      RI1 = IRIG
      SEGACT RI1
C ... NRIG1  = taille du tableau d'infos ...
      NRIG1  = RI1.IRIGEL(/1)
C ... NRIGE1 = nombre de "sous-matrices" ...
      NRIGE1 = RI1.IRIGEL(/2)
      ISYM = 0
      IANT = 0
      IF (NRIG1.GE.7) THEN
         DO 10 IN = 1,NRIGE1
            IANTI = RI1.IRIGEL(7,IN)
            IF (IANTI.EQ.0) ISYM = ISYM + 1
            IF (IANTI.EQ.1) IANT = IANT + 1
            IF (IANTI.EQ.2) THEN
               ISYM = ISYM + 1
               IANT = IANT + 1
            ENDIF
10       CONTINUE
      ELSE
         ISYM = NRIGE1
      ENDIF
*
*     On extrait les matrices ...
*
      NRIGE  = NRIG1
      IF(IRET.EQ.3) THEN
         NRIGEL = ISYM
      ELSE IF(IRET.EQ.4) THEN
         NRIGEL = IANT
      ENDIF
      SEGINI MRIGID
      IRIS = MRIGID
      MTYMAT = RI1.MTYMAT
      ICHOLE = 0
      IMGEO1 = 0
      IMGEO2 = 0
      IFORIG = RI1.IFORIG
      ISUPEQ = 0
      II = 0
      DO 20 IN = 1,NRIGE1
         IF(NRIG1.GE.7) THEN
            IANTI = RI1.IRIGEL(7,IN)
         ELSE
            IANTI = 0
         ENDIF
C    ... Si la matrice possède déjà la bonne symétrie,
C        on ne fait que recopier les pointeurs, sinon,
C        si la matrice est non-symétrique on va en
C        extraire ce qu'il faut ...
         IF (IANTI.EQ.2.OR.IANTI.EQ.IRET-3) THEN
            II = II + 1
            COERIG(II) = RI1.COERIG(IN)
            DO 22 IN2 = 1,NRIGE
               IRIGEL(IN2,II) = RI1.IRIGEL(IN2,IN)
 22         CONTINUE
C       ... si la matrice est non-symétrique, il faut
C           surcharger le flag de symétrie (N° 7) et
C           le pointeur sur IMATRI (N° 4) ...
            IF(IANTI.EQ.2) THEN
               IRIGEL(7,II) = IRET-3
               xMATR1 = RI1.IRIGEL(4,IN)
               SEGACT,xMATR1
               NELRIG = xmatr1.re(/3)
               nligrp = xmatr1.re(/2)
               nligrd= xmatr1.re(/1)
               SEGINI xMATRI
               xmatri.symre=irigel(7,ii)
               IRIGEL(4,II) = xMATRI
               DO 25 INUMEL = 1,NELRIG
*                 XMATR1 = IMATR1.IMATTT(INUMEL)
*                  SEGACT XMATR1
*                  NLIGRD = XMATR1.RE(/1)
*                 NLIGRP = XMATR1.RE(/2)
                  IF(NLIGRD.NE.NLIGRP) THEN
                     write(*,*) 'Matrice non carrée !!!'
                     return
                  ENDIF
*                  SEGINI XMATRI
*                  IMATTT(INUMEL) = XMATRI
                  DO 26 ILIG=1,NLIGRP
                     DO 27 ICOL=1,NLIGRD
                        IF(IRET.EQ.3) RE(ICOL,ILIG,inumel) =
     &    (XMATR1.RE(ICOL,ILIG,inumel)+XMATR1.RE(ILIG,ICOL,inumel))/2.d0
                        IF(IRET.EQ.4) RE(ICOL,ILIG,inumel) =
     &    (XMATR1.RE(ICOL,ILIG,inumel)-XMATR1.RE(ILIG,ICOL,inumel))/2.D0
 27                  CONTINUE
 26               CONTINUE
*                  SEGDES XMATRI
*                  SEGDES XMATR1
 25            CONTINUE
               SEGDES xMATR1,xmatri
            ENDIF
         ENDIF
 20   CONTINUE
 
      SEGDES MRIGID
      SEGDES RI1
      CALL ECROBJ('RIGIDITE',IRIS)
*
      RETURN
      END