C CAPAC3    SOURCE    FANDEUR   11/07/19    21:15:27     7042           

C=======================================================================
C=                            C A P A C 3                              =
C=                            -----------                              =
C=                                                                     =
C=  Fonction :                                                         =
C=  ----------                                                         =
C=  Calcul de la matrice de CAPACITE CALORIFIQUE pour des elements de  =
C=  COQUE TRIANGLE (COQ3) a integration semi-analytique                =
C=                                                                     =
C=  Parametres :  (E)=Entree  (S)=Sortie                               =
C=  ------------                                                       =
C=   NEF      (E)    Numero de l'ELEMENT FINI dans NOMTP               =
C=   IPMAIL   (E)    Numero du segment IMODEL dans le segment MMODEL   =
C=   CLAT     (E)    Chaleur latente du changement de phase            =
C=   IPRIGI   (E/S)  Matrice de CAPACITE (RIGIDITE) resultat (ACTIF)   =
C=======================================================================

      SUBROUTINE CAPAC3 (NEF,IPMAIL,IPINTE,IVAMAT,NVAMAT,IVAPHA,NVAPHA,
     &                   IPMATR,NLIGR)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)

-INC CCOPTIO
-INC CCREEL
-INC CCHAMP

-INC SMCHAML
-INC SMCOORD
-INC SMELEME
-INC SMINTE
-INC SMRIGID

      SEGMENT MMAT1
        REAL*8 XE(3,NBNN),FORME(NBNN)
        REAL*8 CAPSS(NBNN,NBNN),CAPV(NLIGR,NLIGR)
      ENDSEGMENT

      SEGMENT MPTVAL
        INTEGER IPOS(NS),NSOF(NS),IVAL(NCOSOU)
        CHARACTER*16 TYVAL(NCOSOU)
      ENDSEGMENT

C= Quelques constantes numeriques
      PARAMETER (X1s15=0.066666666666666666666666666667D0)
      PARAMETER (X2s15=0.133333333333333333333333333333D0)
      PARAMETER (X8s15=0.533333333333333333333333333333D0)
      PARAMETER (X1s30=0.033333333333333333333333333333D0)

      REAL*8 VACOMP(3)

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 =====
      MPTVAL = IVAMAT
c*      SEGACT,MPTVAL
C- Test sur la constance du champ d'epaisseur : supprime
c*      IPMELV = IVAL(3)
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- Calcul du volume associe a ce point de Gauss (jacobien)
          S1=XZero
          S2=XZero
          S3=XZero
          S4=XZero
          S5=XZero
          S6=XZero
          DO iNoe = 1, NBNN
            S1=S1+SHPTOT(2,iNoe,iGau)*XE(2,iNoe)
            S2=S2+SHPTOT(3,iNoe,iGau)*XE(3,iNoe)
            S3=S3+SHPTOT(3,iNoe,iGau)*XE(2,iNoe)
            S4=S4+SHPTOT(2,iNoe,iGau)*XE(3,iNoe)
            S5=S5+SHPTOT(3,iNoe,iGau)*XE(1,iNoe)
            S6=S6+SHPTOT(2,iNoe,iGau)*XE(1,iNoe)
          ENDDO
          SurfX=S1*S2-S3*S4
          SurfY=S4*S5-S2*S6
          SurfZ=S6*S3-S5*S1
          DJAC = ABS(SurfX*SurfX+SurfY*SurfY+SurfZ*SurfZ)
C- Verification que le volume n'est pas nul en ce point de Gauss
          IF (DJAC.LT.XPETIT) THEN
            INTERR(1) = iElt
            CALL ERREUR(259)
            GOTO 9990
          ENDIF
          DJAC = SQRT(DJAC)
C- Calcul du terme Rho.Cp.Vol en ce point de Gauss
C- (avec prise en compte d'un eventuel changement de phase)
          CALL CAPA5(iGau,iElt,IVAMAT,NVAMAT,IVAPHA,NVAPHA, VACOMP)
          CAPA = DJAC * POIGAU(iGau) * VACOMP(1)
C- Calcul de la contribution du point de Gauss a la matrice
C- CAPACITE elementaire pour cet element fini
          CALL ZERO(CAPSS,NBNN,NBNN)
          do iou=1,nbnn
            forme(iou)=shptot(1,iou,igau)
          enddo
          CALL NTNST(FORME,CAPA,NBNN,1,CAPSS)
C- Ajout de termes specifiques dus a l'integration (analytique)
C- suivant l'epaisseur de l'element de type COQUE
C =======
          EP = VACOMP(3)
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