Sources Esope | Cast3M

tadve8
C TADVE8    SOURCE    CB215821  26/03/06    21:15:09     12485          
 
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,ITYP1)
 
      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
 
-INC TMPTVAL
 
C     Reprise de la MACRO de tadve1.eso (a modifier ensemble) c'est pour ITYP1
      MACRO,(THERMIQUE,DIFFUSION,NAVIER_STOKES,MECANIQUE)
 
      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
 
      SEGMENT,MMAT2
        REAL*8 CEM(NLIGR,NLIGR)
      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
      if (ityp1.eq.MECANIQUE) SEGINI,MMAT2
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)
        if (ityp1.eq.MECANIQUE) CALL ZERO(CEM,NLIGR,NLIGR)
*
* COORDONNEES DES NOEUDS DE L'ELEMENT IEL DANS LE REPERE GLOBAL
*
        CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE)
*
*       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
*
*         Recuperation de GRAD (gradient de FORME et 1er calcul du jacobien)
          CALL TCOND5(IGAU,NBNN,NDIM,XE,SHPTOT,SHP,GRAD,DJAC)
          IF (DJAC.LT.XZERO) IFOIS=IFOIS+1
          DJAC = ABS(DJAC)
          IF (DJAC.LT.XPETIT) THEN
            INTERR(1) = iElt
            CALL ERREUR(259)
            GOTO 9990
          ENDIF
 
*         Recuperation de FORME
          CALL CAPA4(NEF,IGAU,NBNN,XE,SHPTOT,SHP,FORME,DJAC2)
 
          XSCAL = DJAC*POIGAU(IGAU)
 
*         Recuperation des valeurs des composantes (ordre dans tadve1.eso)
*         THERMIQUE & MECANIQUE : XSCAL = DJAC . POIGAU . RHO . C
*         DIFFUSION             : XSCAL = DJAC . POIGAU . CDIF
          DO i = 1, IVAL(/1)
C            IF (IVAL(i).NE.0) THEN
              MELVAL = IVAL(i)
              IBMN = MIN(IEL ,VELCHE(/2))
              IGMN = MIN(IGAU,VELCHE(/1))
 
              if (i.le.IDIM) then
*               Construction du vecteur vitesse (IDIM composante)
                V22(i) = VELCHE(IGMN,IBMN)
 
              else
C               On multiplie les autres scalaires entre eux
                XSCAL = XSCAL * VELCHE(IGMN,IBMN)
              endif
 
C            ELSE
C              if (i.le.IDIM) then
C                V22(i) = XZERO
C              endif
C            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 * XSCAL
          ENDDO
 
 
          IF (ISYMM.EQ.1) THEN
*           CAS SYMETRIQUE : Option 'SYMM' de 'ADVE';
            DO i = 1, NBNN
              r_z = V77(i)
              LI = idim*(i-1)
              DO j = 1, i
                CEL(i,j) = CEL(i,j)
     &                   + (r_z*FORME(j) + V77(j)*FORME(i))/2.D0
                LJ = idim*(j - 1)
                if (ityp1.eq.MECANIQUE) then
                  do k = 1,IDIM
                    CEM(LI + k,LJ + k) = CEL(i,j)
                  enddo
                endif
              ENDDO
            ENDDO
 
          ELSE
*           CAS NON SYMETRIQUE : Par defaut
            DO i = 1, NBNN
              r_z = V77(i)
              LI = idim*(i-1)
              DO j = 1, NBNN
                CEL(j,i) = CEL(j,i) + (r_z *FORME(j))
                LJ = idim*(j - 1)
                if (ityp1.eq.MECANIQUE) then
                  do k = 1,IDIM
                    CEM(LJ + k,LI + k) = CEL(j,i)
                  enddo
                endif
              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
*         CAS SYMETRIQUE
          if (ityp1.eq.MECANIQUE) then
            CALL REMPMT(CEM,NLIGR,RE(1,1,iel))
          else
            CALL REMPMT(CEL,NLIGR,RE(1,1,iel))
          endif
 
        ELSE
*         CAS NON SYMETRIQUE
          if (ityp1.eq.MECANIQUE) then
            CALL REMPMS(CEM,NLIGR,RE(1,1,iel))
          else
            CALL REMPMS(CEL,NLIGR,RE(1,1,iel))
          endif
        ENDIF
      ENDDO
 
C  3 - MENAGE : DESACTIVATION/DESTRUCTION DE SEGMENTS
C ====================================================
 9990   CONTINUE
      SEGSUP,MMAT1
      if (ityp1.eq.MECANIQUE) segsup MMAT2
 
      RETURN
      END