Sources Esope | Cast3M

tadve8
C TADVE8    SOURCE    CB215821  22/04/01    14:03:16     11326          
 
C=======================================================================
C=                            T A D V E 8                              =
C=                            -----------                              =
C=                                                                     =
C=  Fonction :                                                         =
C=  ----------                                                         =
C=  Calcul de la matrice d'ADVECTION pour les                          =
C=  les elements finis MASSIFs a integration NUMERIQUE                 =
C=                                                                     =
C=  Parametres :  (E)=Entree  (S)=Sortie                               =
C=  ------------                                                       =
*     NEF     (E)  NUMERO DE L'ELEMENT-FINI DANS NOMTP (VOIR CCHAMP)
*     IPMAIL   E)  NUMERO DU MAILLAGE ELEMENTAIRE CONSIDERE
*
* AUTEUR, DATE DE CREATION:
* -------------------------
*    MARINO ARROYO, 18 MAI 1999
*
* LANGAGE:
* --------
*     ESOPE + FORTRAN77
*
************************************************************************
      SUBROUTINE TADVE8 (NEF,IPMAIL,IPINTE,IMATE,IVAMAT,NVAMAT,ISYMM,
     &                   IPMATR,NLIGR)
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
 
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
-INC CCHAMP
 
-INC SMCHAML
-INC SMCOORD
-INC SMELEME
-INC SMINTE
-INC SMRIGID
 
      SEGMENT MPTVAL
        INTEGER IPOS(NS),NSOF(NS),IVAL(NCOSOU)
        CHARACTER*16 TYVAL(NCOSOU)
      ENDSEGMENT
 
      SEGMENT,MMAT1
        REAL*8 CEL(NBNN,NBNN),XE(3,NBNN)
        REAL*8 SHP(6,NBNN),GRAD(NDIM,NBNN)
C*      REAL*8 FORME(NBNN),V77(NBNN),V22(IDIM) &lt;- A verifier pour V22
        REAL*8 FORME(NBNN),V77(NBNN),V22(NDIM)
        REAL*8 CMAT(NDIM,NDIM),CMAT1(IDIM,IDIM),CMAT2(IDIM,IDIM)
      ENDSEGMENT
 
C      SEGMENT ,MAXE
C        REAL*8 TXR(IDIM,IDIM),XLOC(3,3),XGLOB(3,3)
C      ENDSEGMENT
 
C  1 - INITIALISATIONS ET VERIFICATIONS
C ======================================
      MELEME = IPMAIL
c*      SEGACT,MELEME
      NBNN   = NUM(/1)
      NBELEM = NUM(/2)
C =====
      MINTE = IPINTE
c*      SEGACT,MINTE
      NBPGAU = POIGAU(/1)
C =====
      MPTVAL = IVAMAT
c*      SEGACT,MPTVAL
C =====
      XMATRI = IPMATR
c*      SEGACT,xMATRI*MOD
c*      NLIGRP = NBNN = NLIGR
c*      NLIGRD = NBNN = NLIGR
C =====
C Recuperation des fonctions de forme et de leurs derivees au
C        centre de gravite de l'element pour le calcul des axes locaux
C        d'orthotropie ou d'anisotropie
C =====
C      IF (IMATE.EQ.2 .OR.IMATE.EQ.3) THEN
C        NLG = NUMGEO(NEF)
C        CALL RESHPT(1,NBNN,NLG,NEF,0,IPINT1,IOK)
Cc*of      IF (IOK.EQ.0) GOTO 999
C        MINTE1 = IPINT1
C        SEGACT,MINTE1
C        NBSH = MINTE1.SHPTOT(/2)
C      ENDIF
 
C =====
C Initialisation des segments de travail
C =====
      IF (IFOMOD.EQ.1) THEN
         NDIM = 3
      ELSE
         NDIM = IDIM
      ENDIF
      SEGINI,MMAT1
C      IF (IMATE.EQ.2 .OR. IMATE.EQ.3) THEN
C        SEGINI,MAXE
C      ENDIF
 
C  2 - BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IPMAIL
C ============================================================
      DO IEL = 1, NBELEM
*
* MISE A ZERO DU TABLEAU CEL
*
        CALL ZERO(CEL,NBNN,NBNN)
*
* COORDONNEES DES NOEUDS DE L'ELEMENT IEL DANS LE REPERE GLOBAL
*
        CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE)
*
CB215821 : En ADVECTION, les vitesses sont donnees dans le repere global
CC Calcul des axes locaux d'orthotropie ou d'anisotropie
C        IF (IMATE.EQ.2 .OR. IMATE.EQ.3) THEN
C          CALL RLOCAL(XE,MINTE1.SHPTOT,NBSH,NBNN,TXR)
C          IF (nbsh.EQ.-1) THEN
C            CALL ERREUR(525)
C            GOTO 9990
C          ENDIF
C        ENDIF
*
*  BOUCLE SUR LES POINTS DE GAUSS
*
        IFOIS = 0
 
        DO IGAU = 1, NBPGAU
*
*           CALCUL DE LA MATRICE GRADIENT DES FONCTIONS DE FORME ET
*           DU JACOBIEN,EN UN POINT DE GAUSS
*
          CALL TCOND5(IGAU,NBNN,NDIM,XE,SHPTOT,SHP,GRAD,DJAC)
          IF (IERR.NE.0) GOTO 9990
          IF (DJAC.LT.XZERO) IFOIS=IFOIS+1
*  Marino calcul de la matrice des fonctions de forme et confirmation du jacobien
          CALL CAPA4(NEF,IGAU,NBNN,XE,SHPTOT,SHP,FORME,DJAC2)
          IF ((ABS(DJAC-DJAC2)).GT.1.d-2) THEN
            WRITE(*,*) '###ERREUR DANS ADVE: Marino jacob diff '
            INTERR(1) = iElt
            CALL ERREUR(259)
            GOTO 9990
          ENDIF
          DJAC = ABS(DJAC)
          IF (DJAC.LT.XPETIT) THEN
            INTERR(1) = iElt
            CALL ERREUR(259)
            GOTO 9990
          ENDIF
          DJAC = DJAC*POIGAU(IGAU)
 
* Recuperation des valeurs des composantes du champ vectoriel
          DO i = 1, IDIM
            IF (IVAL(i).NE.0) THEN
              MELVAL = IVAL(i)
              IBMN = MIN(IEL ,VELCHE(/2))
              IGMN = MIN(IGAU,VELCHE(/1))
              V22(i) = DJAC*VELCHE(IGMN,IBMN)
            ELSE
              V22(i) = XZERO
            ENDIF
          ENDDO
 
C         La vitesse est donnee dans le repere global (elements massifs)
C         Il n'y a pas a distinguer les cas ISOTROPE, ORTHOTROPE et ANISOTROPE
          DO i = 1, NBNN
            r_z = XZERO
            DO j = 1, NDIM
              r_z = r_z + GRAD(j,i)*V22(j)
            ENDDO
            V77(i) = r_z
          ENDDO
 
* CAS SYMETRIQUE
          IF (ISYMM.EQ.1) THEN
            DO i = 1, NBNN
              r_z = V77(i)
              DO j = 1, i
                CEL(i,j) = CEL(i,j)
     &                   + (r_z*FORME(j) + V77(j)*FORME(i))/2.D0
              ENDDO
            ENDDO
* NON SYMETRIQUE
          ELSE
            DO i = 1, NBNN
              r_z = V77(i)
              DO j = 1, NBNN
                CEL(j,i) = CEL(j,i) + (r_z *FORME(j))
              ENDDO
            ENDDO
          ENDIF
        ENDDO
 
C Erreur si, en un point de Gauss, le jacobien change de signe.
        IF (IFOIS.NE.0.AND.IFOIS.NE.NBPGAU) THEN
          INTERR(1) = iElt
          CALL ERREUR(195)
          GOTO 9990
        ENDIF
 
* REMPLISSAGE DE XMATRI
        IF (ISYMM.EQ.1) THEN
          CALL REMPMT(CEL,NLIGR,RE(1,1,iel))
        ELSE
          CALL REMPMS(CEL,NLIGR,RE(1,1,iel))
        ENDIF
C
      ENDDO
*
*     DESACTIVATION DES SEGMENTS
*
C  3 - MENAGE : DESACTIVATION/DESTRUCTION DE SEGMENTS
C ====================================================
 9990   CONTINUE
      SEGSUP,MMAT1
C      IF (IMATE.EQ.2.OR.IMATE.EQ.3) THEN
C        SEGDES,MINTE1
C        SEGSUP,MAXE
C      ENDIF
 
      RETURN
      END