Sources Esope | Cast3M

capac1
C CAPAC1    SOURCE    OF166741  25/02/21    21:15:21     12166          
 
C=======================================================================
C=                            C A P A C 1                              =
C=                            -----------                              =
C=                                                                     =
C=  Fonction :                                                         =
C=  ----------                                                         =
C=  Calcul de la matrice de CAPACITE CALORIFIQUE pour des elements de  =
C=  COQue axisymetrique (COQ2) a integration semi-analytique           =
C=                                                                     =
C=  Parametres :  (E)=Entree  (S)=Sortie                               =
C=  ------------                                                       =
C=   NEF      (E)    Numero de l'ELEMENT FINI dans NOMTP               =
C=   IMAIL    (E)    Numero du segment IMODEL dans le segment MMODEL   =
C=   IPCHA1   (E)    Pointeur sur un segment MCHEL1 de caracteristiques=
C=   CLAT     (E)    Chaleur latente du changement de phase            =
C=   IPRIGI   (E/S)  Matrice de CAPACITE (RIGIDITE) resultat (ACTIF)   =
C=                                                                     =
C=  P. DOWLATYARI, juin 1990 (adaptation de capanu.eso)                =
C=======================================================================
 
      SUBROUTINE CAPAC1 (NEF,IPMAIL,IPINTE,IVAMAT,NVAMAT,IVAPHA,NVAPHA,
     &                   IPMATR,NLIGR,INFOR)
 
      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
 
      SEGMENT MMAT1
        REAL*8 XE(3,NBNN),FORME(NBNN)
        REAL*8 CAPSS(NBNN,NBNN),CAPV(NLIGR,NLIGR)
        REAL*8 VACOMP(NVAMAT)
      ENDSEGMENT
 
C= Quelques constantes numeriques
      PARAMETER (X2Pi=6.283185307179586476925286766559D0)
      PARAMETER (X1s15=0.066666666666666666666666666667D0)
      PARAMETER (X2s15=0.133333333333333333333333333333D0)
      PARAMETER (X8s15=0.533333333333333333333333333333D0)
      PARAMETER (X1s30=0.033333333333333333333333333333D0)
 
      CHARACTER*16 MOFOR
 
 
C- Element purement axisymetrique :
      IF (IFOMOD.NE.0) THEN
        CALL ERREUR(19)
        RETURN
      ENDIF
 
C  1 - INITIALISATIONS ET VERIFICATIONS
C ======================================
      MELEME = IPMAIL
c*      SEGACT,MELEME
      NBNN   = NUM(/1)
      NBELEM = NUM(/2)
      NBNN2 = 2*NBNN
c*      NBNN3 = 3*NBNN
C =====
      MINTE = IPINTE
c*      SEGACT,MINTE
      NBPGAU = POIGAU(/1)
C- Petit test utile ?
      NBNO   = SHPTOT(/2)
      IF (NBNO.NE.2) THEN
        CALL ERREUR(5)
        RETURN
      ENDIF
C =====
      MPTVAL = IVAMAT
c*      SEGACT,MPTVAL
c*C- Verification sur la constance du champ d'epaisseur :
c*C- epaisseur toujours placee en derniere position du mptval
c*      IPMELV = IVAL(NVAMAT)
c*      CALL QUELCH(IPMELV,IOK)
c*      IF (IOK.NE.0) THEN
c*        CALL ERREUR(566)
c*        GOTO 9990
c*      ENDIF
C =====
c*      IF (IVAPHA.NE.0) THEN
c*        MPTVAL = IVAPHA
c*        SEGACT,MPTVAL
c*      ENDIF
C =====
      XMATRI = IPMATR
c*      SEGACT,XMATRI*MOD
c*      NLIGRP = NBNN3 = NLIGR
c*      NLIGRD = NBNN3 = NLIGR
C =====
      SEGINI,MMAT1
 
C  2 - BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IPMAIL
C ============================================================
      DO iElt = 1, NBELEM
C =====
C  2.1 - Recuperation des coordonnees GLOABLES des noeuds de l'element
C =====
        CALL ZERO(CAPV,NLIGR,NLIGR)
C =====
C  2.2 - Recuperation des coordonnees GLOABLES des noeuds de l'element
C =====
        CALL DOXE(XCOOR,IDIM,NBNN,NUM,iElt,XE)
C =====
C  2.3 - Boucle sur les points de Gauss de l'element iElt
C =====
        DO iGau = 1, NBPGAU
C =======
C  2.3.1 - Calcul du volume associe a ce point de Gauss (jacobien)
C =======
          DLX = SHPTOT(2,1,iGau)*XE(1,1)+SHPTOT(2,2,iGau)*XE(1,2)
          DLY = SHPTOT(2,1,iGau)*XE(2,1)+SHPTOT(2,2,iGau)*XE(2,2)
          DJAC = SQRT(DLX*DLX+DLY*DLY)
C- Prise en compte de l'axisymetrie
          CALL DISTRR(XE,SHPTOT(1,1,iGau),NBNN,RR)
          DJAC = ABS(X2Pi*RR*DJAC)
C =======
C  2.3.3 - Verification que le volume (jacobien) n'est pas nul en ce
C          point de Gauss --> Erreur
C =======
          IF (DJAC.LT.XPETIT) THEN
            INTERR(1) = iElt
            CALL ERREUR(259)
            GOTO 9990
          ENDIF
C =======
C  2.3.4 - Calcul du terme Rho.Cp.Vol en ce point de Gauss
C          pour la THERMIQUE
C =======
C          MPTVAL = IVAMAT
          DO i = 1, NVAMAT
            MELVAL = IVAL(i)
            IGMN = MIN(iGau,VELCHE(/1))
            IEMN = MIN(iElt,VELCHE(/2))
            VACOMP(i) = VELCHE(IGMN,IEMN)
          ENDDO
          VALRHO = VACOMP(1)
 
C         CAS THERMIQUE on fait RHO.CP
          IF (INFOR .EQ. 1) VACOMP(1) = VALRHO * VACOMP(2)
 
          CAPA = DJAC * POIGAU(iGau) * VACOMP(1)
C =======
C  2.3.5 - Calcul de la contribution du point de Gauss a la matrice
C          CAPACITE elementaire pour cet element fini
C =======
          CALL ZERO(CAPSS,NBNN,NBNN)
          DO i = 1,NBNN
            FORME(i) = SHPTOT(1,i,iGau)
          ENDDO
          CALL NTNST(FORME,CAPA,NBNN,1,CAPSS)
C =======
C  2.3.6 - Ajout de termes specifiques dus a l'integration (analytique)
C          suivant l'epaisseur de l'element de type COQUE
C =======
          EP = VACOMP(NVAMAT)
C- Erreur si l'epaisseur est est nulle
c*          IF (EP.LE.XPetit) THEN
c*            CALL ERREUR(517)
c*            GOTO 9990
c*          ENDIF
          C1 =  X2s15*EP
          C2 =  X1s15*EP
          C3 = -X1s30*EP
          C4 =  X8s15*EP
          C5 =  C2
          C6 =  C1
          DO j=1,NBNN
            j1 = j + NBNN
            j2 = j + NBNN2
            DO i=1,NBNN
              i1 = i + NBNN
              i2 = i + NBNN2
              Cte = CAPSS(i,j)
              CAPV( i, j) = CAPV( i, j) + C1*Cte
              CAPV(i1, j) = CAPV(i1, j) + C2*Cte
              CAPV(i2, j) = CAPV(i2, j) + C3*Cte
              CAPV(i1,j1) = CAPV(i1,j1) + C4*Cte
              CAPV(i2,j1) = CAPV(i2,j1) + C5*Cte
              CAPV(i2,j2) = CAPV(i2,j2) + C6*Cte
            ENDDO
          ENDDO
        ENDDO
C =====
C  2.4 - Stockage de la matrice de CAPACITE pour cet element fini
C        (remplissage de XMATRI)
C =====
        CALL REMPMT(CAPV,NLIGR,RE(1,1,iElt))
      ENDDO
 
C  3 - MENAGE : DESACTIVATION/DESTRUCTION DE SEGMENTS
C ====================================================
 9990 CONTINUE
      SEGSUP,MMAT1
c*      SEGDES,MELEME,MINTE,XMATRI
c*      MPTVAL = IVAMAT
c*      SEGDES,MPTVAL
c*      IF (IVAPHA.NE.0) THEN
c*        MPTVAL = IVAPHA
c*        SEGDES,MPTVAL
c*      ENDIF
 
      RETURN
      END