C THNUMAC2  SOURCE    PV090527  26/04/30    21:16:41     12529          

C=======================================================================
C=                            T H N U M A C 2                          =
C=                            ---------------                          =
C=                      (TNUMAC dans le cas de la thermique)           =
C=  Fonction :                                                         =
C=  ----------                                                         =
C=  Calcul de la matrice de CAPACITE thermohydrique pour les elements  =
C=  finis MASSIFs a integration NUMERIQUE                              =
C=                                                                     =
C=  Parametres :  (E)=Entree  (S)=Sortie                               =
C=  ------------                                                       =
C=   NEF      (E)   Numero de l'ELEMENT FINI dans NOMTP (cf. CCHAMP)   =
C=   IMAIL    (E)   Numero du segment IMODEL dans le segment MMODEL    =
C=   IPMODE   (E)   Pointeur sur un segment IMODEL suppose ACTIF       =
C=   IPCHEM   (E)   Pointeur sur un segment MCHELM de CARACTERISTIQUES =
C=   IPRIGI  (E/S)  Pointeur sur l'objet RIGIDITE (CONDUCTIVITE)       =
C=                                                                     =
C=   Zakaria HABIBI le 30 juin 2008.                                   =
C=======================================================================

      SUBROUTINE THNUMAC2 (NEF,ipmail,ipinte,ipint1,IVAMAT,NMATT,
     &                     ipmatr,LRE)

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

-INC PPARAM
-INC CCOPTIO
-INC CCREEL

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

-INC TMPTVAL

      SEGMENT MMAT1
        REAL*8 VALMAT(NV1)
        REAL*8 CEL1(NBNN,NBNN),CEL2(NBNN,NBNN),CEL3(NBNN,NBNN)
        REAL*8 CEL4(NBNN,NBNN),CEL5(NBNN,NBNN),CEL6(NBNN,NBNN)
        REAL*8 CEL7(NBNN,NBNN),CEL8(NBNN,NBNN),CEL9(NBNN,NBNN)
        REAL*8 XE(3,NBNN)
        REAL*8 SHP(6,NBNN),GRAD(NDIM,NBNN),FORME(NBNN)
        REAL*8 CMAT(NDIM,NDIM),CMAT1(IDIM,IDIM),CMAT2(IDIM,IDIM)
      ENDSEGMENT

      SEGMENT MAXE
        REAL*8 TXR(IDIM,IDIM),XLOC(3,3),XGLOB(3,3)
      ENDSEGMENT

C INITIALISATIONS ET VERIFICATIONS
C ================================
C Recuperation d'informations sur le maillage elementaire
C =====
      MELEME = IPMAIL
c*      SEGACT,MELEME
      NBNN   = NUM(/1)
      NBELEM = NUM(/2)
C =====
C Recuperation d'informations sur le maillage elementaire
C =====
      MINTE = ipinte
c*      SEGACT,MINTE
      NBPGAU = POIGAU(/1)
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 =====
      IF (ipint1.NE.0) THEN
        MINTE1 = ipint1
c*        SEGACT,MINTE1
        NBSH = MINTE1.SHPTOT(/2)
      ENDIF
C =====
C Initialisation des segments de travail
C =====
      MPTVAL = IVAMAT
      IF (IFOMOD.EQ.1) THEN
        NDIM=3
      ELSE
        NDIM=IDIM
      ENDIF
      NV1 = NMATT
      SEGINI,MMAT1
      MAXE = 0
      IF (ipint1.NE.0) THEN
        SEGINI,MAXE
      ENDIF
C =====
C Matrice de CAPACITE thermohydrique
C =====
      XMATRI = ipmatr
c*      SEGACT,XMATRI*MOD
c*      NLIGRP = 3*NBNN = LRE
c*      NLIGRD = 3*NBNN = LRE

C BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IPMAIL
C ======================================================
      DO iElt=1,NBELEM
C ===
C Mise a zero de la matrice de CAPACITE de l'element iElt
C ===
        CALL ZERO(CEL1,NBNN,NBNN)
        CALL ZERO(CEL2,NBNN,NBNN)
        CALL ZERO(CEL3,NBNN,NBNN)
        CALL ZERO(CEL4,NBNN,NBNN)
        CALL ZERO(CEL5,NBNN,NBNN)
        CALL ZERO(CEL6,NBNN,NBNN)
        CALL ZERO(CEL7,NBNN,NBNN)
        CALL ZERO(CEL8,NBNN,NBNN)
        CALL ZERO(CEL9,NBNN,NBNN)
C ===
C Recuperation des coordonnees GLOBALES des noeuds de l'element
C ===
        CALL DOXE(XCOOR,IDIM,NBNN,NUM,iElt,XE)
C ===
C Calculs des axes locaux d'orthotropie ou d'anisotropie
C ===
        IF (ipint1.NE.0) THEN
          CALL RLOCAL(XE,MINTE1.SHPTOT,NBSH,NBNN,TXR)
          IF (NBSH.EQ.-1) THEN
            CALL ERREUR(525)
            GOTO 9990
          ENDIF
        ENDIF
C ===
C Boucle sur les points de Gauss de l'element iElt
C ===
        iFois=0
        DO iGau=1,NbPGau
C - Calcul du jacobien, des fonctions de forme et de leurs
C   derivees au point de Gauss iGau
          CALL CAPA4(NEF,iGau,NBNN,XE,SHPTOT,SHP,FORME,DJAC)
          DJAC=ABS(DJAC)
          IF (IERR.NE.0) GOTO 9990
          IF (DJAC.LT.XZero) iFois=iFois+1
C - Erreur si le jacobien est nul en ce point de Gauss
          IF (DJAC.LT.XPetit) THEN
            INTERR(1)=iElt
            CALL ERREUR(259)
            GOTO 9990
          ENDIF
          DJAC=DJAC*POIGAU(iGau)
C - Recuperation de la ou des valeurs de conductibilite au point
C   de Gauss iGau (tableau VALMAT)
          DO i=1,NMATT
            IF (IVAL(i).NE.0) THEN
              MELVAL=IVAL(i)
              IBMN=MIN(iElt,VELCHE(/2))
              IGMN=MIN(iGau,VELCHE(/1))
              VALMAT(i)=VELCHE(IGMN,IBMN)
            ELSE
              VALMAT(i)=XZERO
            ENDIF
          ENDDO
C - Cas d'un materiau ISOTROPE de conductibilite K
C   Calcul de la contribution du point de Gauss a la matrice
C   CAPACITE elementaire de cet element fini
C   Ajout du terme XK*transposee(N)*N
C   Seul cas valide en dimension 1
             XK=VALMAT(10)*DJAC
             CALL NTNST(FORME,XK,NBNN,1,CEL1)
             XK=VALMAT(11)*DJAC
             CALL NTNST(FORME,XK,NBNN,1,CEL2)
             XK=VALMAT(12)*DJAC
             CALL NTNST(FORME,XK,NBNN,1,CEL3)
             XK=VALMAT(13)*DJAC
             CALL NTNST(FORME,XK,NBNN,1,CEL4)
             XK=VALMAT(14)*DJAC
             CALL NTNST(FORME,XK,NBNN,1,CEL5)
             XK=VALMAT(15)*DJAC
             CALL NTNST(FORME,XK,NBNN,1,CEL6)
             XK=VALMAT(16)*DJAC
             CALL NTNST(FORME,XK,NBNN,1,CEL7)
             XK=VALMAT(17)*DJAC
             CALL NTNST(FORME,XK,NBNN,1,CEL8)
             XK=VALMAT(18)*DJAC
             CALL NTNST(FORME,XK,NBNN,1,CEL9)

        ENDDO
C =====
C Erreur si, en un point de Gauss, le jacobien change de signe
C =====
        IF (iFois.NE.0.AND.iFois.NE.NbPGau) THEN
          INTERR(1)=iElt
          CALL ERREUR(195)
          GOTO 9990
        ENDIF
C =====
C Stockage de la matrice de CAPACITE pour cet element fini
C =====  (remplissage de XMATRI)
       CALL REMPMCH
     & (CEL1,CEL2,CEL3,CEL4,CEL5,CEL6,CEL7,CEL8,CEL9,NBNN,RE(1,1,ielt))

      ENDDO

C  3 - MENAGE : DESACTIVATION/DESTRUCTION DE SEGMENTS
C ====================================================
 9990 CONTINUE
      SEGSUP,MMAT1
      IF (ipint1.GT.0) THEN
        SEGSUP,MAXE
      ENDIF

      RETURN
      END